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Abstract 

Heat  transfer  with  change  of  phase  (freezing  or  melting)  is  important  in 
numerous  scientific  and  engineering  applications.  Since  the  pioneering  works 
of  Lam6  and  Clapeyron,  Neumann  and  Stefan,  o  number  of  anolytical  and 
numerical  techniques  have  been  developed  to  deal  with  freezing  end  melting 
problems.  One  such  analytical  tool  is  the  method  of  perturbation  expansions, 
which  is  the  main  focus  of  this  work.  The  report  begins  with  a  review  of  the 
perturbation  theory  and  outlines  the  regular  perkirbation  method,  the  method  of 
strained  coordinates,  the  method  of  matched  asymptotic  expansions,  and  the 
recently  developed  method  of  extended  perturbation  series.  Next,  the  applications 
of  these  techniques  to  phase  change  problems  in  Cartesian,  cylindrical,  and 
spherical  systems  are  discussed  in  detail.  Although  the  bulk  of  the  discussion 
is  confined  to  one-dimensional  situations,  the  report  also  includes  two-  and 
three-dimensional  cases  where  admittedly  the  success  of  these  techniques  has 
so  far  been  limited.  The  presentation  is  sufficiently  detailed  that  even  the  reader 
who  is  unfamiliar  with  the  perturbation  theory  can  understand  the  material. 
However,  at  the  same  time,  the  discussion  covers  the  latest  literoture  on  the 
subject  and  therefore  should  serve  as  a  state-of-the-art  review. 
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NOMENCLATURE 


a  aiiq>litude  of  temperature  oscillation,  °C  (°F)  or  a  constant 

Oj  j  binomial  coefficients 

ag  series  coefficients 

A  amplitude  of  temperature  oscillation,  dimensionless  or  area  of  cross  section, 

m^  (ft^). 

Ao  constant 

A,„  series  coefficients 

b  a  constant 

Bq  a  function  of  Bi  and 

bi  a  function  of  Bi  and  y  or  a  constant 

63  a  function  of  rf 

b„  series  coefficients 

B|  a  constant 

Bi  convective  Biot  number 

Bi'r  radiative  Biot  number 

c  specific  heat,  kJ/kg  K  (Btu/lbm®F)  or  a  constant 

c  average  specific  heat,  kJ/kg  K  (Btu/lbm°F) 

Cf  specific  heat  at  temperature  Tf,  kJ/kg  K  (Btu/lbm°F) 

Cn  series  coefficients 

C,Ci,C2,C3  constants 

e  surface  emissivity  for  radiation,  dimensionless 

e\  first  Shanks  transformation 

E  freezing  front  location,  dimensionless 

/  function  of  e 

fiB)  function  of  0 

n,  functions  of  Xf 

F  temperature  or  overall  radiation  shape  factor,  dimensionless 
F„  temperature  series  coefficients 

g  velocity  of  freezing  front,  dimensionless 

Si  coefficients  of  series  for  g 

h  convective  heat  transfer  coefficient,  W/(m^  K)  (Btu/ft^hr°F)  or  an  integer 
ij  integers 

k,k^  thermal  conductivity  of  solid  phase,  W/m  K  (Btu/ft  hr°F) 

kf  thermal  conductivity  at  temperature  Tf 

ki  thermal  conductivity  of  liquid  phase,  W/m  K  (Btu/ft  hr°F) 
k  average  thermal  conductivity,  W/m  K  (Btu/ft  hr°F) 

K  ratio  k^kf  or  a  constant 

Ki  a  constant 

t  slab  thickness,  m  (ft) 

L  latent  heat  of  fusion,  kJ/kg  (Btu/lbm) 

Li  a  constant 

m,n  integers 

P  perimeter  of  cross  section,  m  (ft) 

q  heat  flux,  W/m^  (Btu/ft^hr) 

Q  heat  flux,  dimensionless 

r  radial  coordinate,  dimensionless 

rf  radial  location  of  freezing  front,  dimensionless 
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R  radial  coordinate,  m  (ft) 

Rf  radial  location  of  freezing  front,  m  (ft) 

R^  wall  radius,  m  (ft) 

s  Laplace  variable 

5  Stefan  number,  dimensionless 

Vl.-^n^n+l  partial  sums 

t  time,  sec 

T  temperature,  “C  (®F) 

Tg  ambient  temperature,  °C  ("F) 

Tf  freezing  temperature,  “C  (°F) 

Ti  initial  temperature,  ®C  (®F) 

Tq  coolant  temperature  or  mean  coolant  temperature  or  wall  temperature  or  fin  base 
temperature.  °CC*F) 

Tv  vaporization  temperature,  “C  (“F) 

u  temperature  distribution  in  the  solid  phase,  dimensionless 
Ui,  u„  coefficients  of  series  for  u 

V  temperature  distribution  in  the  liquid  phase  or  steady-state  velocity  of  vaporizing 

front,  dimensionless 
v„  coefficients  of  series  for  v 

X  distance  from  wall,  m  (ft) 

X*  stretched  distance,  m  (ft) 

Xf  location  of  freezing  front,  m  (ft) 

Xg  reference  of  characteristic  distance,  m  (ft) 

Xy  location  of  vaporizing  front,  m  (ft) 

X  distance  from  wall,  dimensionless 

Xf  location  of  freezing  front,  dimensionless 

Xfh  coefficients  of  series  for  Xf 

z  time,  dimensionless 

Greek 

a  thermal  diffusivity,  m^/s  (ft^/s)  or  slope  of  the  specific  heat-temperature  curve 
divided  by  Cf,  (°F~*)  or  ambient  temperature,  dimensionless  or  slope  of 
Domb-Sykes  plot  or  a  =  Vaj/ttf 
Oj  thermal  diffusivity  of  solid  phase,  m^/s  (ft^/s) 
a/  thermal  diffusivity  of  liquid  phase,  m^/s  (ft^/s) 

P  inverse  of  Stefan  number  5  or  the  ratio  of  convective  and  radiative  Biot  numbers  or 

the  slope  of  thermal  conductivity-temperature  curve  divided  by  kf,  (®C"’  or 
oF-i) 

E,£i  ,82  perturbation  parameters  or  coordinates 
Eq  radius  of  convergence  of  series 

e*  Eulerized  perturbation  parameter 

q  distance  or  similarity  variable  or  inverse  of  rf  or  vaporization  velocity,  all  dimen¬ 
sionless 

q„  coefficients  of  series  for  q 

q*  q  for  inner  expansion 

qH  coefficients  of  series  for  q* 

qf  similarity  variable  based  on  the  thermal  diffusifity  of  liquid  phase 

q,  similarity  variable  based  on  the  thermal  diffusivity  of  solid  phase 

6  temperature,  dimensionless 
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6a  ambient  teperature,  dimensionless 

6f  freezing  temprature,  dimensionless 

6i  initial  temperature,  dimensionless 

6/  tenq)erature  of  the  liquid  phase,  dimensionless 
6a  temperature  of  the  solid  phase,  dimensionless 
6*  6  for  the  inner  expansion 

6o  coefficients  for  series  for  6* 

X  freezing  or  vaporizing  front  location,  dimensionless 

eigenvalues 

coefficients  of  series  for  X 

X^  location  of  stationary  freezing  front,  dimensionless 

X*  X  for  inner  expansion 

p.  fin  parameter,  m“'  (fr*) 

p  density,  kg/n?  (Ibm/ft^) 

a  Stefan-Boltzmann  constant  or  location  of  vaporization  front,  dimensionless 
Gi  straining  functions 

a„  coefficients  of  series  for  a  (the  location  of  vaporization  front) 

a*  inner  expansion  variable  for  o  (die  location  of  vaporization  front) 

oJI  coefficient  of  series  o* 

X  time,  dimensionless 

Xp  preheating  time,  dimensionless 

%  coefficients  of  series  for  x 

Xy  time  to  complete  vaporization,  dimensionless 

X*  X  for  inner  expansion  or  x  after  Euler  transformation 
frequency  of  temperature  oscillation,  sec~^ 

<0  frequency  of  temperature  oscillation,  dimensionless 
<|>  strained  variable,  or  temperature  (outer  expansion),  dimensionless 
<|>*  temperature  (inner  expansion),  dimensionless 

<t»n  coefficients  of  series  for  (ji* 

i|r  strained  variable 


Perturbation  Techniques  in  Conduction-Controlled 
Freeze-Thaw  Heat  Transfer 

VIRGIL  J.  LUNARDINI  AND  ABDUL  AZIZ 


1.  INTRODUCTION 

Heat  transfer  with  freezing  or  melting  occurs  in  a  number  of  applications  such  as  ice  formation, 
permafrost  melting,  metal  casting,  food  preservation,  storage  of  latent  energy,  automatic  welding,  etc. 
The  vast  literature  that  exists  on  the  subject  has  been  codified  recently  into  a  monograph  by  Lunardini 
( 1 99 1 ).  In  the  simplest  models,  the  two  phases  in  a  typical  freezing  or  melting  problem  are  separated 
either  by  a  sharp  boundary  or  a  mushy  zone  that  moves  with  time.  Since  the  location  of  this  moving 
boundary  or  zone  is  not  known  a  priori,  the  solution  is  difficult  to  achieve  even  when  the  heat  transfer 
process  is  assumed  to  be  conduction  controlled.  Most  exact  solutions  rely  on  similarity  transforma¬ 
tions,  introduced  by  Lamd  and  Clapeyron  ( 1 83 1 ),  or  Neumann  ( 1 860),  and  have  been  examined  in 
some  detail  (Lunardini  1 98 1 , 1 99 1 ).  It  is  only  under  highly  idealized  conditions  that  one  can  obtain 
useful,  exact,  analytical  solutions. 

For  more  realistic  situations,  one  must  therefore  think  in  terms  of  either  an  approximate  analytical 
technique  or  a  fully  numerical  approach  based  on  finite  differences  or  finite  elements.  The  former 
strategy  is  particularly  convenient  for  preliminary  calculations.  Over  the  years,  several  approximate 
techniques  have  been  devised  to  solve  freezing  and  melting  problems,  including  perturbation  methods, 
the  heat  balance  integral  technique,  and  variational  methods. 

This  report  is  devoted  to  the  application  of  perturiiation  techniques  to  freezing  and  melting 
problems.  The  coverage  begins  with  an  overview  of  the  perturbation  theory  and  includes  a  brief 
discussion  of  the  basic  concepts,  regular  perturbation  method,  singular  perturbation  techniques  such 
as  the  method  of  strained  coordinates  and  the  method  of  matched  asymptotic  expansions,  and  the 
method  of  extended  perturbation  series.  The  applications  of  these  methods  will  cover  both  one¬ 
dimensional  and  two-dimensional  phase  change  problems  in  Cartesian,  cylindrical  and  spherical 
systems.  The  accuracy  of  the  perturbation  solutions  will  be  assessed  in  the  light  of  other  approximate 
solutions  and  also  of  exact  and  numerical  results. 

Asymptotic  methods  are  excellent  for  the  study  of  freeze-thaw  problems.  Singular  perturbation 
techniques  may  be  applied  to  treat  the  singularities  that  occur  when  phases  appear  or  disappear.  In 
principle,  multidimensional  problems  can  be  handled  and  perturbation  methods  are  often  superior  for 
dealing  with  nonlinear  boundary  conditions  or  convective  effects.  Perturbation  also  yields  valuable 
insights  into  the  basic  physics  of  the  problems.  The  main  disadvantage  of  the  method  is  the  increasing 
difficulty  of  obtaining  higher  order  terms. 


2.  REVIEW  OF  PERTURBATION  THEORY 

The  method  of  perturbation  expansion  is  a  well  established  analytical  tool  that  has  found 
applications  in  many  areas  of  engineering.  The  subject  is  covered  in  detail  in  several  currently 


available  books.  Those  dealing  with  engineering  applications  include  Nayfeh  (1981),  Van  Dyke 
( 1 975)  and  Aziz  and  Na  ( 1 984).  While  the  first  two  discuss  applications  in  solid  mechanics  and  fluid 
dynamics,  the  1n"t  one  is  exclusively  devoted  to  problems  in  heat  transfer. 

In  freezing  and  melting  problems,  the  main  difficulty  is  the  presence  of  a  moving  boundary  that 
separates  the  solid  and  liquid  phases,  but  several  other  difficulties  can  arise.  These  additional 
difficulties  may  be  due  to  nonlinear  or  time-dependent  boundary  conditions,  finite  phase-change 
domain,  an  intermediate  mushy  zone,  domain  geometry,  etc.  Some  of  these  difficulties  can  be  handled 
by  the  perturbation  mi’‘'tod,  as  will  be  demonstrated  later  in  the  report. 

Perturbation  jry  is  based  on  the  concept  of  an  asymptotic  solution.  If  the  basic  equations 
describing  a  phase-change  problem  can  be  expressed  such  that  one  of  the  parameters  or  variables  is 
small  (or  very  large)  then  the  full  equations  can  be  approximated  by  letting  the  perturbation  quantity 
approach  its  limit  and  an  approximate  solution  can  be  found  in  terms  of  this  perturbation  quantity.  Such 
a  solution  approaches  a  limit  as  the  perturbation  quantity  approaches  zero  (or  infinity)  and  is  thus  an 
asymptotic  solution.  The  result  can  often  be  improved  by  expanding  in  a  series  of  successive 
approximations,  the  first  term  of  which  is  the  limiting  solution.  One  then  has  an  asymptotic  series  or 
expansion.  Thus  we  perturb  the  limiting  solution  by  parameter  or  coordinates. 

One  is  then  concerned  with  the  asymptotic  expansion,  generally  for  a  small  parameter  such  as  the 
Stefan  number,  of  the  solutions  of  the  conduction  equation  with  solidification. 

The  first  step  in  a  perturbation  analysis  is  to  identify  the  perturbation  quantity.  This  is  done  by 
expressing  the  mathematical  model  in  a  dimensionless  form,  assessing  the  order  of  magnitude  of 
different  terms  and  identifying  the  term  that  is  small  compared  to  others.  The  coefficient  of  this  term, 
which  could  be  a  dimensionless  parameter  or  a  dimensionless  variable,  is  then  chosen  as  a  perturbation 
quantity  and  designated  by  the  symbol  e.  Once  e  is  identified,  the  solution  is  assumed  as  an  asymptotic 
series  of  e.  Next,  this  series  solution  is  substituted  into  the  governing  equations  for  the  problem.  By 
equating  the  coefficients  of  each  power  of  e  to  zero,  one  can  generate  a  sequence  of  subproblems.  These 
problems  are  solved  in  succession  to  obtain  the  unknown  coefficients  of  the  series  solution. 

The  foregoing  procedure  is  termed  parameter  perturbation  or  coordinate  perturbation  depending 
on  whethere  is  aparameteror  acoordinate.  In  either  case,  a  further  distinction  is  made  between  regular 
perturbation  if  the  expansion  is  uniformly  valid  and  singular  perturbation  if  the  expansion  fails  in 
certain  regions  of  the  domain.  When  a  singular  perturbation  expansion  is  encountered,  the  usefulness 
of  the  solution  is  limited  unless  it  can  be  rendered  uniformly  valid.  Note  that  the  terms  in  the  expansion 
need  not  be  convergent  for  the  results  to  be  useful  since  its  asymptotic  nature  assures  that  only  a  few 
terms  may  yield  adequate  accuracy  for  small  values  of  E. 

The  two  main  techniques  for  achieving  uniform  validity  that  have  been  used  in  freezing  and  melting 
problems  are  the  method  of  strained  coordinates  and  the  method  of  matched  asymptotic  expansions. 
In  the  method  of  strained  coordinates,  both  the  dependent  and  the  independent  variables  are  expanded 
in  terms  of  e  such  that  the  coefficients  of  the  two  series  are  functions  of  new,  unknown  independent 
variables.  The  assumed  series  expansions  are  substituted  into  the  governing  equations,  and  the 
unknown  coefficients  are  found  by  ensuring  that  higher  approximations  are  no  more  singular  than  the 
first  one.  The  procedure  leads  to  an  implicit  but  uniformly  valid  solution. 

The  method  of  matched  asymptotic  expansions  achieves  uniform  validity  by  supplementing  the 
regular  perturbation  expansion,  which  is  now  called  the  outer  expansion,  with  an  inner  expansion  in 
which  the  independent  variable  is  stretched  out  such  that  it  describes  the  behavior  in  the  region  where 
the  outer  expansion  breaks  down.  A  uniformly  valid  solution  is  finally  derived  by  matching  the  outer 
and  the  inner  expansion  according  to  Van  Dyke’s  or  some  other  matching  principle. 

In  most  instances  the  perturbation  expansions  are  terminated  at  the  second  or  the  third  term.  There 
are  two  reasons  for  such  an  early  truncation.  First,  the  values  of  e  of  interest  are  often  small  compared 
to  unity,  and  therefore  the  truncated  expansion  that  converges  rapidly  for  small  values  of  e  is 
sufficiently  accurate.  The  second  reason  is  that  the  higher-order  terms  of  the  series  are  increasingly 
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more  difficult  to  calculate.  Sometimes,  however,  it  may  be  desirable  to  improve  the  solution  so  that 
it  can  be  used  for  values  of  e  that  are  not  too  small.  This  can  be  achieved  with  the  method  of  extended 
perturbation  series.  The  approach  consists  of  three  steps.  The  first  step  is  to  write  a  computer  program 
to  solve  the  sequence  of  perturbation  equations  either  in  symbolic  form  or  numerically  and  generate 
a  large  number  of  terms.  The  second  step  is  to  examine  the  coefficients  generated  in  the  first  step  to 
identify  the  singularities  in  the  complex  plane  that  restrict  the  convergence  of  the  series.  This 
information  about  the  analytic  structure  of  the  solution  is  finally  used  to  transform  the  series  using  one 
or  a  combination  of  techniques  such  as  the  Euler  transformation,  Pad6  approximants.  Shanks 
transformation,  series  reversion,  etc.  The  three-step  procedure  results  in  a  series  solution  that  is 
accurate  over  a  wider  range  of  e  than  the  original  series. 


3.  REGULAR  PERTURBATION  EXPANSIONS 

Ttie  regular  perturbation  method  has  proved  to  bean  effective  tool  to  solve  a  variety  of  freezing  and 
melting.  Most  applications  involve  parameter  perturbations  but  as  shown  later,  the  coordinate 
perturbation  approach  is  also  feasible  in  some  cases.  The  method  will  be  illustrated  with  the  help  of 
a  number  of  examples. 

3.1  The  Stefan  problem 

The  mathematical  class  of  problems  for  which 
a  moving  boundary  exists  are  often  called  Stefan 
problems.  The  name  derives  from  Stefan’s  ( 1891) 
early  work  on  sea  ice  formation.  In  this  section 
the  Stefan  problem  will  be  limited  to  the  follow¬ 
ing  case.  The  classical  Stefan  problem  considers 
the  freezing  of  a  saturated  liquid  of  semi-infinite 
extent  as  shown  in  Figure  1 .  The  liquid  is  as¬ 
sumed  to  be  initially  at  its  freezing  temperature 
Tf.  At  the  time  t  ~  0"^,  the  surface  temperature  at 
X  =  0  is  suddenly  reduced  to  a  subffeezing  value 
Tg  <  Tf.  The  lowering  of  the  surface  temperature 
causes  the  liquid  to  freeze.  With  time,  the  freez¬ 
ing  front  progresses  in  the  direction  of  increas¬ 
ing  X  values.  Let  .Xf  be  the  thickness  of  the  solid 
phase  at  any  instant  of  time.  The  unfrozen  liquid 
continues  to  remain  at  Tf  throughout  the  solidification  process.  With  appropriate  thermal  property 
values,  the  Stefan  problem  can  also  describe  the  melting  of  a  solid  of  semi-infinite  extent  initially  at 
its  melting  temperature  and  heated  by  a  higher  wall  temperature. 

The  Stefan  problem  is  a  highly  idealized  model  for  the  actual  freezing  or  melting  process.  The 
mathematical  model  to  be  described  is  based  on  the  following  assumptions: 

1 .  Both  solid  and  liquid  phases  are  homogeneous  and  isotropic. 

2.  The  phase  change  occurs  at  a  discrete  temperature  and  consequently  there  is  no  mushy  zone 
containing  a  mixture  of  two  phases. 

3.  The  two  phases  are  separated  by  a  sharply  defined  interface  (front)  instead  of  whisker-like 
(dendritic)  growth  observed  experimentally  (Sparrow  et  al.  1979). 

4.  Conduction  is  the  only  heat  transfer  mode. 

For  the  case  of  freezing,  the  temperature  distribution  in  the  solid  phase  is  uescribed  by  the  following 
mathematical  model; 


Figure  /.  One-dimensional  freezing  of  a  semi¬ 
infinite  region. 
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(1) 


d^T 

r(0./)=7’o.  T{xf,t)=T{  (2) 


k 


a;c 


X=Xf 


=  pL^ 
dt 


(3) 


where  k,  p  and  a  represent  the  thermal  conductivity,  density,  and  thermal  diffusi  vity,  respectively  of 
the  solid  phase.  The  quantity  L  is  the  latent  heat  of  solidification. 

To  prepare  tlie  foregoing  model  for  a  perturbation  analysis,  we  recognize  that  the  Stefan  number, 
which  signifies  the  importance  of  sensible  heat  to  the  latent  heat,  is  small  in  some  phase  change 
applications  and  therefore  it  can  serve  as  a  perturbation  quantity.  For  example,  water/ice  systems  or 
soil  systems  typically  have  Stefan  numbers  less  than  1/2,  unless  the  boundary  temperatures  are  very 
high.  With  this  in  mind,  we  introduce  the  following  dimensionless  quantities: 

T  =  kt/pcx}  ,£=c(T{-  tJ/L  (4) 


where  c  is  the  specific  heat  of  the  solid  phase  and  is  a  reference  distance.  Using  eq  4  to 
nondimensionalize  eq  1, 2, 3  we  obtain 


_ae 

32,2  "ax 


(5) 


0(o.x)=l,  0U  =  Xf,T)=O,  — 

ax 


t  dx 


(6) 


The  perturbation  quantity  e  which  appears  in  the  boundary  condition  at  the  interface  can  be  brought 
into  the  governing  equation  (eq  5)  by  noting  that  the  transient  storage  term,  a0/ax,  will  be  related  to 
the  transient  motion  of  the  interface,  dX^dx.  Thus  we  can  think  of  0  as  a  function  of  two  new 
independent  variables,  X,  the  nondimensional  distance  as  before,  andXf,  the  instantaneous  location  of 
solid/liquid  interface.  For  the  time  derivative,  then,  we  obtain 


30  _  30  aXf 
3x  3Xf  dz 

From  eq  6  we  have 


(7) 


dXj 

dt  dx 


X=X( 


Using  eq  8  in  eq  7  and  the  result  in  eq  5,  the  transient  heat  conduction  equation  becomes 


3  0  =  ^ 


3X' 


3Xf 


30 


3X 


X  =  Xf 


The  first  two  boundary  conditions  (eq  6)  can  now  be  written  in  terms  of  X  and  Xf  as 


(8) 


(9) 
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e(x=o,  Xf)  =  i,  e{x  =  Xf.  Xf)=o  (lO) 

Let  us  derive  a  perturbation  solution  by  assuming  an  asymptotic  series  solution  of  the  form 

e  =  i  0nt”  (11) 

n  =  0 

Retaining  only  the  first  three-terms  of  the  series  (eq  11),  we  have 

9  =  0Q  +  £0 1  +  £2®2  (12) 

Substituting  eq  12  into  eq  9  and  eq  10,  we  get 

^^00  +  +  £2  ^^02  _ 

3x2  ^^2 

_e[i§o +gi^  +  e2^1  +e^  +  (13) 

[3Xr  3Xr  3X J  [  3X  3X  3X 

eo(X=0,Xf)+  £01  (X=0,Xf)+  £202(x=O,Xf)  =  1  (14) 

0o(X  =  Xf,Xf)+ £10,  (X=Xf,  Xf)+ £202(X  =  Xf,Xf)  =  O  (15) 


0o(X  =  Xf,Xf)+ £10,  (X=Xf,  Xf)+ £202(X  =  Xf,Xf)  =  O  (15) 

By  equating  coefficients  of  like  powers  of  £  in  eq  13,  Hand  15  we  obtain  the  following  sub-problems. 
Zero-order,  e° 


o 

1) 

rS  ^ 

(16) 

9o{X=0,Xf)=l, 

0o(x  =  Xr,Xf)=O 

(17) 

First-order, 

3^0 I  _  300  300 

3X2  "  axf  3X 

X=Xf 

(18) 

0l(X=O,  Xf)  =  0, 

0i(X  =  Xf.Xf)  =  O 

(19) 

Second-order,  ^ 

d\  _  [300  301 
3X2  [axf  ax 

^  300  30| 

X=X(  X=X(^^^. 

(20) 

02(X=O,  Xf)  =  0, 

02(X  =  Xf,  Xf)  =0 

(21) 
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It  may  be  noted  that  the  zeroth-order  problem  represents  the  quasi-steady  approximation  to  the 
problem  which  was  used  by  Stefan  (Lunardini  1981). 

Integrating  eq  16  twice  and  applying  the  boundary  conditions  (eq  17),  the  solution  for  Bq  is  found 
to  be 

00=1-  i  (22) 

Using  eq  22  on  the  right-hand  side  of  eq  1 8  gives 


=JL 
dX^  x} 

whose  solution,  subject  to  eq  19,  is 


(23) 


(24) 


Now  we  use  eq  22  and  eq  24  on  the  right  hand-side  of  eq  20  to  obtain 

dX^  [6  \xn  2  /. 

which  can  be  integrated  twice  using  eq  21  to  give 


(25) 


02  = 


— W 

360  \X(I 


Using  eq  22, 24,  and  26  in  eq  1 2,  the  three-term  perturbation  solution  takes  the  form 


(26) 


(27) 


90 

The  equation  for  Xf  can  now  be  derived  by  evaluating  - — 

oX 

result  on  the  right-hand  side  of  eq  8.  The  resulting  equation  is 


from  eq  27  and  substituting  the 


X=Xf 


=  J-  ( e  -  J-  +  X  e3j  +  0  (e'*) 

A  Xf  '  3  45  ' 

Integrating  eq  28  subject  to  the  initial  condition,  t  =  0,  Xf  =  0,  we  get 

X^  =  2t  (e-  +  -2-  e^j 

1  3  45  / 

If  we  wish  to  express  T  as  a  function  of  Xf,  we  write  eq  29  as 

X  =  J- X^  £■' ( 1  - -L  £  + -^  e^]  ' 

2  \  3  45  f 


(28) 


(29) 


(30a) 
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Using  binomial  expansion 


fl-ie+-Le^)  *  =  l+i£--2-e^+0  (e^) 
\  3  45  I  3  45 


leads  to 


t  =  J-e"*  e — L  e  x^+0(e^) 

2  6  45 


(30b) 


By  using  the  similarity  technique,  the  exact  analytical  solution  of  eq  5,  6  can  be  obtained  as 
[Lunardini  1991,Ozisik  1980] 


erf(XX/Xf) 

er^ 

where  X  =  Xf/2Vt  is  the  root  of  the  equation 
Vn  Xe^  erfX  =  e 


(31) 


(32) 


It  is  interesting  to  demonstrate  that  the  perturbation  solution  eq  29  reproduces  the  exact  solution  (eq 
32)  for  sufficiently  small  values  of  X.  Consider  the  series  expansions  for  e^  and  erf  X  valid  for  small 
values  of  X: 


e^^  =  1  +  ^  +  . . . 

2  6 


erfX  =  -2-[x-  ^  +  ...] 

Vtt  \  3  10  42  / 

Substituting  eq  33  and  eq  34  into  eq  32  and  simplifying,  we  obtain 

(2  x2+4  x4+ J_x6+-l^X*+  ...)=e 
V  3  15  105  / 


(33) 

(34) 


(35) 


A  general  series  2  ^  ®  can  be  reversed  and  written  as 

n  s  1 


X2=  X  *n£" 

n=  1 


(36) 


where  the  first  three  coefficients  b„,  for  example,  are 


*1 


=  -L.t2  =  -«2.i3  =  a£2Zfif3 


a\ 


(37) 


Noting  from  eq  35  that  a\  =  2, 02  =  4/3, 03  =  8/15,  we  find  from  eq  37  that  b\  =  1/2, 62  =  -1/6.  and  63 
=  7/90.  Using  these  values  in  eq  36,  we  find 
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=  -L  e  -  -L  (38) 

2  6  90 

Now  =  x}lAx  and  hence  eq  38  can  be 
written  as 

g2  ^  7  g3 

4x  2  6  90 

or 

X^=2tfe- le2  +  J-e3)  (39) 

V  3  45  / 

which  matches  exactly  with  the  perturbation 
solution  eq  29. 

In  Figure  2,  the  three-term  perturbation 
solution  for  the  freezing  time  t,  eq  30  is 
compared  with  the  exact  analytical  solution 
eq  32  for  e  =  0.2, 0.4  and  0.6.  The  agreement 
is  excellent  even  at  e  =  0.6. 

3..2.  Planar  freezing  of  a  saturated  liquid 
with  convective  cooling 

The  problem  to  be  considered  here  is  a 
variation  of  the  Stefan  problem  of  the  pre¬ 
vious  section.  The  constant  temperature 
boundary  condition  at  x  =  0  is  now  replaced 
with  a  convective  boundary  condition  as 
shown  in  Figure  3.  The  coolant  temperature  is 
Tg  and  die  convective  heat  transfer  coefficient 
is  h.  The  problem  is  described  by  eq  1,  2,  3 
except  that  the  first  boundary  condition  in  eq 
2  now  becomes 


Figure  2.  Perturbation  solution  for  Stefan  problem. 


Figure  3.  Freezing  of  semi-infinite  saturated 
medium  with  surface  convection. 


k^{0,t)=h[T(0,t)-T^  (40) 

dx 

In  this  case,  it  is  more  convenient  to  introduce  the  following  dimensionless  quantities: 


T(-  To  k  pck  L 


The  resulting  dimensionless  model  is 


dX^  ~  dXf  dX 


^(0,Xf)=e{0,Xf)-l 


0(X  =  Xf,  Xf)=0. 


(43) 
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Assuming  a  three-tenn  solution  of  the  form  of  eq  12,  and  following  the  procedure  of  section  3.1 , 
we  find  that  the  governing  equations  for  6o,  0].  and  62  are  still  given  by  eq  16-21  except  that  die  wall 


boundary  conditions  are  now  given  by 

^(aXf)=0o(O.Af)-l  (44) 

dX 

^(O.Xf)=0,(O,Xf)  (45) 

oX 

^(O,Xf)=02(O.Xf)  (46) 

aX 

The  solution  of  eq  16  subject  to  eq  44  and  eq  17  is 

0o  =  — (^f-X)  (47) 

1  +  Xf 

Similarly,  the  solution  of  eq  1 8  subject  to  eq  45  and  eq  19b  is 

0 1  = - 1 - [  ( 1  +  Xf )  X^  +  3  ( 1  +  Xf )  X^  -  ( 3  +  Xf )  x/  X  -  ( 3  +  Xf )  Xf^]  (48) 

£  f  1  . 


Finally  the  solution  for  02  subject  to  eq  46  and  eq  21b  can  be  obtained  as 

0  2  = - 1 - [  ( 1 9  xf  +  1 1 4  X?  +  225  Xf  +  360  Xf )  ( 1  +  X) 

360(l+Xf)’ 

- 10  Xf  (1  +  Xf)  (Xf  +  3  Xf  +  12)  ( 3  +  X)X^ 

-9  (1+ Xf)2  (5+ X)X'*]  (49) 

Following  the  procedure  of  section  3. 1 ,  the  solution  for  freezing  time  t  can  be  obtained  as 

t=i  e-i[(l  +  Xf)2-l]+  -_L^[{l+Xf)3-3(l  +  Xf)  +  2] 

2  6(1  +  Xf) 

-  J-e - ! - [(l+Xf)®-5(l+Xf)^  +  9(l+Xf)-5]  (50) 

45  (l+Xf)^ 

The  perturbation  analysis  presented  in  this  section  was  first  discussed  by  Pedroso  and  Domoto 
(1973a)  who  presented  a  scheme  for  programming  the  perturbation  calculations  so  that  a  digital 
computer  can  be  used  to  generate  as  many  terms  as  desirable.  The  scheme  circumvents  the  mounting 
algebraic  labor  entailed  in  manual  calculations  as  one  goes  to  higher  order  terms.  They  present  their 
results  for  t  as 


ee 

eT=  X 

fi  =  0 

where  the  coefficients  Cn  are  functions  of  Xf. 


(51) 
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Table  1.  Numerical  values  of  c,. 


Xff 

1 

2 

3 

4 

5 

6 

7 

8 

9 

0.2 

0.2200 

0.01778 

-41.001564 

0.0003945 

-0.0001383 

0.00005842 

-0.00002793 

0.00001462 

-0.0000)8213 

0.4 

0.4800 

0.06476 

-0.008154 

0.002846 

-0.001335 

0.0007332 

-0.0004438 

0.0002869 

-0.0001940 

0.6 

0.7800 

0.1350 

-0.01932 

0.00748 

-0.003804 

0.002211 

-0.001389 

0.0009138 

0.0006180 

0.8 

1.120 

0.2252 

-0.03398 

0.01368 

-0.007083 

0.004134 

-0.002573 

0.001661 

-0.001096 

1 

1.500 

0.3333 

-0.05139 

0.02092 

-0.01083 

0.006253 

-0.003828 

0.002422 

-0.001576 

1.4 

2.380 

OA980 

-0.09282 

0.02768 

-0.01918 

0.01082 

-0.006461 

0.004001 

-0.002643 

1.8 

3.420 

0.9257 

-0,1418 

0.05692 

-0.02852 

0.01585 

-0.009353 

0.005750 

-0.004005 

2.2 

4.620 

1.311 

-0.1979 

0.07854 

-0.03893 

0.02145 

-0.01259 

0.007719 

-0.005770 

2.6 

5.980 

1.753 

-0.2608 

0.1025 

-0.05048 

0.02767 

-0.01620 

0.009923 

-0.007997 

3.0 

7.500 

2.250 

-0.3305 

0.1292 

-0.06323 

0.03456 

-0.02021 

0.01237 

-0.01071 

4 

12.00 

3.733 

-0.5347 

0.2068 

-0.1006 

0.05479 

-0.03199 

0.01956 

-0.01963 

5 

17.50 

5.555 

-0.7823 

0.3008 

-0.1459 

0.07939 

-0.04632 

0.02832 

-0.03153 

Table  1  gives  the  first  nine  values  of  c„  for  a  range  of  values  of  Xf.  The  first  three  values  agree  with 
those  given  by  eq  50. 

The  present  problem  was  also  solved  by  Huang  and  Shih  ( 1 975)  by  first  immobilizing  the  interface 
position  using  the  Landau  transformation  and  then  using  a  regular  perturbation  analysis.  Their  three- 
term  solutions  for  the  temperature  distribution  and  freezing  time  match  exactly  with  the  present  results. 
The  use  of  the  Landau  transformation  makes  the  nonlinearity  due  to  the  moving  interface  explicit  and 
facilitates  the  subsequent  use  of  the  perturbation  method.  However,  the  approach  adopted  here 
achieves  the  same  simplicity  without  a  formal  introduction  of  the  Landau  transformation. 

The  present  problem  is  a  special  case  of  a  more  general  problem  treated  by  Westphal  (1967)  who 
considered  the  freezing  of  a  semi-infinite  medium  cooled  by  the  coolant  having  a  time-dependent 
temperature,  T(y{t).  Westphal’ s  exact  solution  is  in  the  form  of  integrals  and  infinite  series,  which  make 
it  extremely  awkward  for  numerical  computations.  Another  exact  solution  has  been  presented  by 
Lozano  and  Reemsten  ( 1 98 1 ),  but  again  the  solution  is  tedious  if  numerical  information  is  to  be  derived 
from  it.  In  view  of  the  complexity  of  these  exact  solutions,  the  perturbation  solution  is  best  compared 
with  the  heat  balance  integral  solution  of  Goodman  (1958),  analog  solution  of  Krieth  and  Romie 
(1955),  and  the  local  similarity  solution  of  Aziz  and  Lunardini  (1991).  The  comparison  shown  in 


□  Perturbation 
O  HBI  (Goodman.  1958) 

A  Analog  (KreitharKf  Romie,  1955) 


Figure  4.  Approximate  solution  for 
freezing  with  surface  convection. 


S  =  cAT/L 
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Figure  4  demonstrates  that  the  prediction  of  the  perturbation  solution  is  consistent  with  die  other 
approximate  solutions. 


33  Planar  freezing  of  a  saturated  liquid  due  to 
sinusoidal  surface  temperature  variation 

This  problem  is  of  considerable  importance  in  modeling  the  cyclical  behavior  of  the  active  layer 
in  permafrost  regions  because  the  actual  annual  ground  surface  temperature  variation  often  approxi¬ 
mates  a  sine  curve  (Kane  et  al.  1 99 1 ,  Aziz  and  Lunardini  1992).  The  problem  is  a  variation  of  the  Stefan 
problem;  the  mathematical  description  is  given  by  eq  1-3  except  that  the  constant  temperature 
boundary  condition  T(0,  t)  =  Tq  is  now  replaced  by  the  following  sinusoidal  condition 

T(0,t)  =  T^  +  as.mQl<Tf  (52) 

Introducing  the  following  dimensionless  quantities, 

0  =(7-ro)/(rf-ro),  X  =  jc(Q/a)''^.  x  =£2r  (53) 

A  =a/(rf-To),  t=c{Tf-T^)lL 
into  eq  1, 2, 3, 52  leads  to 


8^9  99  90 

BX(  dX 


(54) 


0(X  =0,  Xf)  =  Asinx,  0(X  =  Xf,Xf)  =  l  (55) 

Assuming  a  perturbation  expansion  in  the  form  of  eq  12  and  considering  the  first  two  terms  only,  we 
find  that  is  governed  by  eq  1 6  while  0i  is  governed  by  eq  1 8  without  the  negative  sign  on  the  right- 
hand  side.  However,  the  boundary  conditions  change  to 


0o(X=O,  Xf)  =  A  sin  X,  0o(X  =  Xf,  Xf)  =1 
0l{X=O,  Xf)  =  0,  0i(X  =  Xf,  Xf)=  0 
The  solution  of  eq  16  subject  to  eq  56  is 


(56) 

(57) 


00  =  ^  +  ( 1  -  A  sin  X 

Xf  1  Xfl 

Similarly,  the  solution  of  eq  18  subject  to  eq  57  is 


0i  =  1  ( 1  -  A  sin  x)^  1  - 

6  Xf  \Xf  j  . 

Thus  the  two-term  perturbation  solution  becomes 


0=^+1 

(,_X) 

1  A  sin  X  +  J-  e  ( 1  -  A  sin  x)^  X 

[l-(. 

Xj^l 

Xf) 

'  6  Xf 

ff/  J 

(58) 


(59) 


(60) 
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Using  eq  60  to  evaluate,  substituting  it  in  eq  6,  and  noting  that  the  negative  sign  on  the  right-hand 
side  is  to  be  ignored,  we  obtain  the  differential  equation  for  Xf  as  follows: 


dXf 

dt 


=  e 


1  -  A  sin  T 
Xf 


£  ( 1  -  A  sin  t) 
3 


(61) 


Integrating  eq  61  subject  to  the  initial  condition,  t  =  0,  Xf  =  0,  we  obtain 


X,?  =  2e{(x  +  Acost-a)-  |(x  +  2  A  cos  x-2A)  +  ^  A^|t  -  i  sin  2Tjj|  (62) 

If  A  =  0,  the  present  problem  reduces  to  the  Stefan  problem  discussed  in  section  3.1.  The 
perturbation  approach  adopted  for  this  problem  was  first  introduced  by  Lock  et  al.  (1969)  and  Lock 
(1969,  1971),  and  they  reported  that  the  agreement  between  the  perturbation  and  finite  difference 
solutions  was  excellent. 


3.4.  Variable-property  Stefan  problem 

The  variable-property  Stefan  problem  in  which  the  thermal  conductivity  and  specific  heat  of  the 
solid  phase  are  temperature-dependent  has  been  studied  by  Aziz  (1978)  and  Pedroso  and  Domoto 
(1973b)  using  a  regular  perturbation  method.  While  Pedroso  and  Domoto  use  the  Stefan  number  5  as 
a  perturbation  parameter,  Aziz  identifies  a  new  perturbation  parameter  that  is  a  measure  of  the 
variation  of  a  thermal  property  with  temperature  which  makes  the  solution  valid  for  all  values  of  S. 
Both  approaches  will  be  illustrated  here. 

3.4. ].  Pedroso  and  Domoto 's  solution 

Consider  a  semi-infinite  region  of  liquid  initially  at  its  freezing  temperature  Tf.  At  time  r  >  0  the 
face  at  A  =  0  is  maintained  at  constant  temperature  <  Tf.  The  analysis  takes  into  consideration  the 
property  variation  in  the  solid  phase,  the  unfrozen  liquid  being  assumed  to  remain  at  the  freezing 
temperature.  The  governing  equations  are 


pc(r)|:  = 
dt 


dx 


(63) 


7’(0,r)=7’o. 


T{xf,t)=Tf, 


x=jrf 


dt 


By  introducing  the  following  quantities: 


(64) 


k{T)dt 

jTq 

k(T)dt 

T-f-To 

Tf-To 

f  k(T)dt 

JTc 

,  _  c  {Tf- To) 

k{Tf-  To) ’ 

C 

L 
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(65) 


x=^Jllzlo)i,  /(e)=£l^.  T,  =  x 
pL  k{Q)/k  ■'f 

it  can  be  shown  that  eq  63, 64  reduce  to  the  following  boundaiy-value  problem: 


dr\  dt\ 


(66) 


e(0)  =  0,  0(l)l  =  l,  x}  =  2x^ 

d[\ 


(67) 


n=i 


We  solve  eq  66  and  67  assuming /(d)  =  a  +  bQ,  where  a  and  b  are  constants.  For  a  three-term 
expansion  of  the  form  0  =  0o  +  £0i  +  e^02,  the  governing  equations  for  d,,,  dj  and  02  are 


gO.  d\-Q 

df\^ 

(68) 

CD 

o 

O 

It 

o 

o® 

II 

(6ga) 

£> :  +  q(a  +  b0o)  ^  ^ 

(69) 

dr^  dr\  <t(\ 

n=» 

o 

II 

o 

II 

o 

(69a) 

ddi 

d6i 

dr\^ 

<h\ 

dn 

n=i  ^ 

dq 

n=i. 

(o+  bdo) 


+  bQ,^  ^ 


02(O)  =  O,  02(1)  =  0. 

The  above  sequence  of  equations  can  be  integrated  successively  to  give 
00  =  ^ 

01  =  i  Tij^|a+  ^  -  ccr^-  ^  bVj 

02  = - 1 —  ( 2394a2  +  765  +  2111  ab)r\ 

45,360 


(70) 

(71) 


(72) 

(73) 


+  -J-  a{a+  b)n3  +  -1—  b^n^  +  -L  aZfiS 
36  144  40 


+  -1-  abq®  +  -5-  b^q’^ 
30  504 


(74) 
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Thus  the  complete  three-term  perturbation  solution  is 


0  =  11  +  i  eri 
6 


—  p2 


J_ 


45,360 


[(a+i6)-mi2-i6Ti3] 

{2394a2  +  765  *2  +  2772  ab)i\ 


+  -L  a(a+  b)r\^  +  -1—  b^r\^  +  -L  nS 
36  144  40 


+  abx\^  +  -5—  ^2  ^7 
30  504  J  • 

Utilizing  eq  75  in  eq  67,  the  freezing  front  location  is  given  by 


1-  e(^  a  +  2.  b] 

e2(  774 

a2  + J45. 

^2  +  21X 

\3  2  1 

\2835 

3024 

630  n 

(75) 


(76) 


The  accuracy  of  eq  75  and  76  has  been  checked  against  the  direct  numerical  solution  of  eq  66  and 
67  by  Pedroso  and  Domoto  (1973b)  and  found  to  be  good.  Furthermore,  these  authors  give  the 
solutions  for  any  arbitrary  form  of  the  functiony(0). 


3.4.2  Aziz’s  Solution 

Aziz  (1978)  considered  the  solution  of  eq  63  and  64  for  two  cases  of  thermal  property  variation; 
1 )  a  linear  specific  heat-temperature  variation  with  constant  thermal  conductivity;  2)  a  linear  thermal 
conductivity-temperature  variation  with  constant  specific  heat.  Mathematically  these  variations  are 


expressed  as 

c  =C{[l+  a{T(-T)]  (77) 

ifc  =  *f[l+ p(rf-7')]  (78) 

Introducing  eq  77  and  78  into  eq  63  and  64  and  using  the  well  known  similarity  transformation,  the 
governing  equations  in  dimensionless  form  become 

F"  +  2n(l  +  eF)F' =  0  (79) 

for  the  temperature-dependent  specific  heat  and 

(l  +  eF)F‘'+  e(F'f  +  2TlF'=  0  (80) 


for  the  temperature-dependent  thermal  conductivity.  The  boundary  conditions  common  to  eq  79  and 
80  are 


F(0)=1,  F(A,)=0,  F'(X)=-  2XIS 
where 


(81) 


F  =  {Tf-  T)/{Tf-  To),  11  =-^,  X  =  .  ttf  =  ^ ,  S  =  a  (rf-  To),  e  =  afTf-  To) 

2ya{t  2iatt  L 
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for  eq  79,  e  =  p(Tf-  T^)  for  eq  80.  TTie  primes  denote  differentiation  with  respect  to  q. 


Since  e  is  small  for  most  applications  (Aziz  and  Benzies  1976),  F  and  "k  are  expanded  in  the  form 


F  =  i  e"  F„ 


«= 0 


(82) 


X  =  X  (83) 

n=0 

Substituting  eq  82  and  83  into  eq  79  and  80  and  equating  coefficients  of  like  powers  of  e,  a  set  of 
equations  for  F,,,  F]  etc.  can  be  generated.  The  condition  F(0)  =  1  in  eq  81  becomes 

f='o(0)  =  l  Fi(0)=0  /=l,2,..../i  (84) 

However,  the  freezing  front  conditions  in  eq  8 1  contain  £  implicitly.  To  remove  the  implicimess,  we 
expand  F„(X)  and  F  n(X)  in  a  Taylor’s  series  about  to  recast  these  in  explicit  form  as 


X  £" 

n  =0 


I 
(  =  0 


I  e^  Xj 
;=i 


dn' 


(85) 


.  u 

T  £'>  y  d''*'F„(K) 

n  =  0  (=0  l!  dq' 


i  e"Xn  . 


(86) 


Since,  as  shown  later,  the  two-term  perturbation  expansion  gives  accurate  results  for  a  wide  range  of 
parameters,  we  formulate  only  the  zero-  and  first-order  problems  and  present  their  solutions. 
Temperature-dependent  specific  heat.  The  zero-order  problem,  which  is  the  classical  Neumann 


problem,  has  the  solution 

Fo  =  1  -  erfq/erfXo 

(87) 

and 

Vn  Xo  exp(Xo)  erfXo  =  S 

(88) 

The  first-order  problem  is 


F,+2qF,'=-2qFoF, 


(89) 


F,(0)  =  0 


F,(Xo)=-X,Fo(Xo) 


Fo{Xo)  =  -X,F;(Xo)-2X,/S 


(90) 


Using  the  method  of  variation  of  parameters,  the  solution  of  eq  89  and  90  is  finally  obtained  as 
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^  X.o{exp(2Xo)-ll  +  25XoX.i 

crfXo 

-ilexpMf|-rfj^j - 1 - j,_„p(.2,2)) 

^erfXo  \  erfXj  jtlerfA.,;}^ 

X  =Mii  +  1  -  exp  (2XS)} 

'  S{s+2  (1+ S)A.2} 


(91) 


(92) 


Temperature-dependent  thermal  conductivity.  The  solution  of  the  zero-order  problem  is  given  by 
eq  87  and  88.  The  first-order  problem  is 

F,+  2tiF;  =  -FoFo- (fof  (93) 

with  boundary  conditions  given  by  eq  90.  The  method  of  variation  of  parameters  gives  the  solution 
as 


Xi(l-exp(2X^))  +  2SXoX.t+ 
_ 2 

5  ^  erf  Xq 


erf„+5f5PH!!) 

■fit  erf  Xo 


etfT)  ] 

erfXo/ 


and 


+ 


_ ! _ 

rt(erfXo)^ 


_  Xo[5^-  25Xi+25XS{exp(2X^)-l  }] 
‘  25{S-K2(l  +  5)Xi) 


(94) 


(95) 


Aziz  evaluated  the  above  perturbation  series  for  5  =  0.0822, 0.3564, 0.9205, 1 .9956, 4.0601  and 
8. 1720  (Xo  =  0.2, 0.4, 0.6, 0.8, 1.0, 1 .2,  respectively)  with  e  ranging  in  each  case  from -3  to  -h3. 

First  we  discuss  the  results  for  a  temperature-dependent  specific  heat.  For- 1  ^e<  1  the  perturbation 
solutions  agree  to  within  0.5%  with  the  numerical  solutions.  Indeed,  the  convergence  is  so  rapid  that 
even  when  £  is  as  large  as  ±3  the  error  does  not  exceed  1.7%.  A  sample  perturbation  temperature 
distribution  is  given  in  Table  2.  The  results  for  the  freezing  front  location  appear  in  Table  3.  These 


Table  2.  Temperature  distribution  for 5  =  4.0601  (X,  s  1):  effect  of  a 
temperature^lependent  specific  heat. 


n 

c=  1.0 

OJ 

0 

-OJ 

-1.0 

0 

1.0000 

t.OOOO 

1.0000 

1.0000 

1.0000 

0.2 

0.7038 

0.7198 

0.7357 

0.7517 

0.7677 

0.4 

0.4436 

0.4676 

0.4916 

0.5157 

0.5397 

0.6 

0.2374 

0.2604 

0.2834 

0.3064 

0.3295 

0.8 

0.0858 

0.1026 

0.1194 

0.1362 

0.1529 

0.9595 

0.9798 

1 

1.0202 

1.0404 

0 

0 

0 

0 

0 
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Table  3.  Interface  location  parameter  X  at  different  S:  effect  of  a 
temperature-dependent  specific  heat. 


e 

S  =  0.0822 

0.3564 

0.9205 

1.9956 

4.0601 

8.1720 

-1.0 

0.1993 

0.3953 

0.5853 

0.7735 

0.9595 

1.1472 

-0.5 

0,1997 

0.3976 

0.5932 

0.7868 

0.9798 

1.1736 

-0.0 

0.2000 

0.4000 

0.6000 

0.8000 

1.0000 

1.2000 

-0.5 

0.2003 

0.4023 

0.6068 

0.8132 

1.0202 

1.2264 

-1.0 

0.2006 

0.4047 

0.6136 

0.8264 

1.0404 

1.2528 

results  indicate  that  at  low  values  of  5  the  effect  of  a 
variation  of  specific  heat  on  temperature  and  freez¬ 
ing  front  location  is  small  but  becomes  progressively 
significant  as  S  increases.  The  same  conclusion  holds 
for  other  types  of  boundary  conditions,  such  as  a  wail 
heat  flux  varying  linearly  or  exponentially  with  time 
(Chung  and  Yeh  1976)  or  the  wall  being  cooled  by 
combined  convection  and  radiation  (Yeh  and  Chung 
1977). 

In  the  case  of  a  temperature-dependent  thermal 
conductivity,  the  temperature  series  agree  to  within 
2.2%  with  the  numerical  solutions  in  the  range  -0.5 
<  e  <  0.5.  The  maximum  error  of  2.2%  occurs  in  the 
neighborhood  of  the  interface  for  5  =  8. 1 720.  Although  not  reported  here,  a  second-order  correction 
extends  the  range  of  validity  of  the  solution  to  higher  values  of  e.  Fortunately,  the  more  important 
freezing  front  series  converges  rapidly  and  with  first-order  correction  gives  results  within  1 .8%  in  the 
range  -2  ^  e  ^  2.  Sample  results  for  temperature  distribution  and  freezing  front  location  are  given  in 
Tables  4  and  5,  respectively .  Compared  with  a  variable  specific  heat,  the  effect  of  a  variation  in  thermal 
conductivity  is  much  more  pronounced  on  both  the  temperature  and  the  freezing  front  location  at  all 
values  of  S.  This  relatively  strong  effect  of  variable  thermal  conductivity  is  also  evident  in  the  results 
of  Chung  and  Yeh  ( 1 976)  who  analyzed  the  case  of  a  time-dependent  (linear  and  exponential)  wall  heat 
flux  using  a  heat  balance  integral.  It  can  also  be  seen  in  another  paper  (Yeh  and  Chung  1 977)  in  which 
they  study  the  case  of  combined  convective-radiative  cooling  of  the  wall  using  Biot’s  variational 
principle. 


Table  4.  Temperature  distribution  for 
S  =  4.0601  (Xg  =  1):  effect  of  a  temperature- 
dependent  thermal  conductivity. 


n 

e  =  0J 

0 

-OJ 

0 

1.0000 

1.0000 

1.0000 

0.2 

0.8096 

0.7357 

0.6618 

0,4 

0.5960 

0.4916 

0.3872 

0.6 

0.3824 

0.2834 

0.1843 

0.8 

0.1935 

0.1194 

0.0453 

0.9082 

0 

1.0 

0 

1.0918 

0 

Table  5.  Interface  location  parameter  X  at  different  5:  effect  of  a 
temperature-dependent  thermal  conductivity. 


€ 

S  =  0.0822 

0.3564 

0,9205 

1.9956 

4.0601 

8.1720 

1.0 

0.2493 

0.4949 

0.7335 

0.9634 

1.1836 

1.3946 

0.5 

0.2247 

0.4474 

0.6668 

0.8817 

1.0918 

1.2973 

0 

0.2000 

0.4000 

0.6000 

0.8000 

1.0000 

1.2000 

0.5 

0.1753 

0.3526 

0.5332 

0.7183 

0.9082 

1.1027 

1.0 

0.1507 

0.3051 

0.4665 

0.6366 

0.8164 

1.0054 

3.5.  Planar  freezing  with  sinusoidal  variation  of  coolant  temperature 

This  problem  is  also  a  variation  of  the  Stefan  problem  of  section  3. 1 .  The  constant  temperature 
boundary  condition  at  x  =  0  is  now  replaced  by  a  boundary  condition  characterizing  the  sinusoidal 
variation  of  coolant  temperature  which  can  be  expressed  as 
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=  ^0+  asin£2/ 


(96) 


where  is  the  coolant  temperature  fluctuating  around  a  mean  value  of  Tg,  with  an  amplitude  a  and 
frequency  Q.  The  convective  boundary  condition  is  written  as 

*  ^  (0,/)  =  //[r(0,/)-ra(/)]  (97) 

dx 

Equation  3  and  the  second  condition  of  eq  2  also  apply  here. 

The  problem  will  be  solved  using  a  double  series  expansion  involving  two  perturbation  parameters. 
Cl  and  62-  As  defined  below,  ej  represents  the  classical  Stefan  number,  while  £2  is  a  measure  of  coolant 
temperature  fluctuations.  The  problem  is  recast  into  dimensionless  form  by  defining  the  following 
quantities: 


0={7’f-r)/(7’f-7’o),  X=f!^,  ei  =  c{Tf-T^)/L 

k 

£2  =  To),  <0  =  pLkdlh'^iTf-  To)  (98) 

x  =  h^{Tf-To)  tIpLk 
which  give 


a^e  30 

'  ax 


(99) 


ax 


=  0(0, x)  -  1  +  62  sin  ®x 

x=o 


(100) 


0(Xf,x)  =  O 


(101) 


dXf  _  ao 
dz  ax 


X  =  Xt 


(102) 


To  solve  eq  99-102,  a  double  series  expansion  for  0  and  Xf  is  assumed.  Retaining  only  the  first 
three  terms  of  the  series,  we  have 


0  =  0o(X,x)+  ei0i  (X,x)+  e202(X.x) 


(103) 


X{o(x)+  €|Xf,{x)+  e2Xf2(x)  (104) 

Substituting  eq  103  and  104  into  eq  99-102  and  removing  the  implicitness  of  Cj  and  £2  Xf  by 
Taylor  series  expansions  about  X^,  the  following  system  of  equations  for  Oq,  0i ,  02.  Xf,,,  Xfi  and  Xf2 
is  obtained: 


^  =  0  (105) 

ax  2 
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=  eo(0.x)  - 1 ;  0o(Xfo,  T))  =0;  ^ 

at  oX 


8^01  _  ^ 
3X2  -  3^ 


(106) 


(107) 


301 

3X 


Jf=0 


=  0,  (0,t);  0,  (Xfo,  X))  =  -Xfi  QoiXfo,  x); 


dXfi  _  _  30) 

~  w 


X=Xfo 


-  Xfi^ 

3X 


X=Xt 


3^02  ^ 
3X2 


302 

3X 


=  02(0, x)  +sin  cox;  02 (^fo.  t)  =  -Xf2  0o(^fo.  t); 


x=o 


dXf2  _  302 

dt  ~~  dX 


-X(2 


X=Xto 


3X 


X=Xf 


where  primes  denote  differentiation  with  respect  toX. 

The  solutions  of  eq  105-1 10  are  found  to  be  as  follows: 


(108) 


(109) 


(110) 


00  =  (Xfo -X)  Xfo  (111) 

Xfo  =(l+2x)‘^-l  (112) 


+  lXfoXfoX2- 
2 


(113) 


Xf,  =- 


Xi^oXfo 

2(1  + Xfo) 


02  = 


X  +  1-  (l+2x) 


1/2 


(l+2x)'^ 


sin  (OX  - 


(l+y) 

a)(l+2x)02 


( 1-cos  (ox) 


(114) 


(115) 


Xf2  =  -(  1-cos  (ox)/o)( I+2x)'^  (116) 

The  solutions  given  by  eq  1 1 1-1 16  stem  from  the  work  of  Gutman  ( 1 986)  and  have  been  correct¬ 
ed  to  account  for  the  second  term  in  both  of  the  equations  for  Xf  ■  and  Xf  2-  These  terms  were  missed 
by  Gutman. 
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3.6.  Freezing  of  a  semi-infinite  region  with 
convective  and  radiative  cooling 

Another  variation  of  the  Stefan  problem  arises  when  the  boundary  condition  at  jc  =  0  is  that  of 
simultaneous  convection  and  radiation.  A  particular  case  was  considered  by  Yan  and  Huang  (1979) 
for  a  finite  phase  change  region  when  the  other  end  was  insulated.  Since  the  presence  of  radiation 
introduces  a  nonlinear  boundary  condition,  they  linearized  it  by  using  the  approximation 
r  '*  =  4Ta^  T -3  where  is  the  common  ambient  temperature  for  both  convection  and  radiation. 
With  this  approximation,  the  problem  was  reduced  to  that  of  pure  convective  cooling,  the  case 
discussed  in  section  3.2.  Yan  and  Huang  developed  three-term  perturbation  solutions  for  the 
temperature  distribution  and  the  freezing  front  location  and  found  that  the  perturbation  solutions 
agreed  quite  well  with  variational  solutions. 

Here  we  approach  the  semi-infinite  region  case  using  the  double  series  expansion  method  of  the 
previous  section.  Equations  I  -3  apply  except  that  the  first  condition  in  eq  2  is  replaced  by  aeon  vecti  ve- 
radiative  condition 

A:=0,  k^  =  h{T-T^)+aF{T'^-  T^)  (117) 

dx 

To  normalize  the  equations,  we  introduce  the  following  dimensionless  quantities: 
e=T/{T(-T^),  X  =  x/x„  x  =  k{Tf-T^)t/pLxi, 


e\  =  c{T{  -T^IL ,  e2  =  aF(rf-raP  x^lk  ,Bi  =  hx^lk 
into  eq  1-3  and  1 17  and  obtain 

^  =e  ^ 

~  ‘8Yf  ■  dX  „  „ 


(118) 


(119) 


X  =  0,  —  =  fii(e-0a)  +  e2(0'‘-ea) 
dx 

X  =  Xf,  0=0f 
dYf  _  30  I 


(120) 

(121) 


Assuming  an  expansion  for  0  in  the  form  of  eq  103  and  carrying  out  the  procedure  outlined  there, 
we  obtain 

^  =0  (122) 
3X2 

X  =  0,  ^  =  Bi(0o-0a);  X  =  Xf,  00  =  0,  (123) 


3X 


320,  _  300  300 

3X2  axf  3X 


(124) 
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(125) 


^=i?i0,;  X=Xf,  01=0 


^=0 

9X2 

x  =  0,  ^  =  0102+  0o-0o:X=Xf.  02  =  0 
9X 

The  solutions  for  0^,  Oj  and  02  are  as  follows: 


0o  =  0a+(0f-0a)(-^±^) 

\Xf+Bi~'l 

0,  =  i  ■  ( IjLSlX  ]  (xf  +  30/"*  Xf )  -  (x^  +  30/~'  X^) 

1  Xf+0i'*  J  f\t+BiXfJ 


Utilizing  eq  1 28-1 30  to  evaluate  in  eq  1 2 1 ,  the  differential  equation  for  Xf  is  obtained  as 

^  ^  I  e,  [  __0« —  (Xf3+  30,-’  xf)  -  ( 3xf  +  60r’  Xf) 

Xf  +  BC^  6  +  L 1 +BfXf 

+  e2([0a-^(-^^~^‘‘^^*-T-el|  —i—  (131) 

IL  Xf+0r‘  J  / 

For  specified  values  of  0f_  Og,  Bi,  e,,  and  62,  eq  131  can  be  integrated  numerically  with  the  initial 
condition  x  =  0,  Xf = 0  to  give  Xf  as  a  function  of  x.  Figure  5  shows  the  variation  of  the  freezing  front 


—  Equation  (131) 

0.8-  D  Chung  and  Yeh  (1975) 
o  Yan  and  Huang  (1979) 


e, .  4.0 


£1.1.0 


e,  -  3.0 
Bi.1.0 


X'-X/£ 

Figure  5.  Perturbation  solution  for  freeing  with  surface  convection  and  radiation. 
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XfWith  t'=  x/ei  =  atlxl  for0f=4,0a  =  3,5t=  I, £2  =  0.016 and ej  =  0.025, 0.25  and  1.0.  These 
values  were  chosen  so  that  the  present  solution  can  be  compared  with  the  variational  results  of  Chung 
and  Yeh  (1975)  and  a  different  perturbation  solution  reported  by  Yan  and  Huang  (1979).  The 
predictions  of  the  present  perturbation  solution  are  quite  close  to  the  other  reported  solutions. 


4.  REGULAR  COORDINATE  PERTURBATION  EXPANSIONS 

Compared  with  the  parameter  perturbation  expansion,  the  use  of  a  coordinate  perturbation 
approach  to  freezing  and  melting  problem  has  been  very  limited.  Two  examples  will  be  considered, 
one  that  has  been  developed  by  the  present  authors  and  is  being  reported  for  the  first  time,  and  a  second 
one  that  has  been  presented  by  Rathjen  and  Jiji  (1970). 

4.1.  Planar  freezing  of  a  saturated  liquid  with  convective  cooling 

We  revisit  the  problem  of  section  3.2  and  develop  a  coordinate  perturbation  solution  valid  for  short 
times  when  convection  is  strong  or  for  long  times  when  convection  is  weak.  To  solve  the  problem 
which  is  defined  by  eq  1,  3  and  40  and  the  second  condition  of  eq  2,  we  introduce  the  similarity 
transformation  as  follows; 


F=(rf-r)/(rf-ro).  ti  =  Jc/(2VS7)  (i32) 

and  define  a  dimensionless  coordinate  e  =  klilhioit).  It  can  be  shown  that  the  similarity  transfor¬ 
mation  (eq  132)  reduces  eq  1  to  an  ordinary  differential  equation  in  Fas 

F"+2qF'=0.  (133) 

The  boundary  conditions  at  x-0  and  x=Xf  reduce  to 

tl  =  0,  eF  =  F-l;  Ti  =  A.  =  Xf/(2Va7),  F=0; 

q=X,  F'=-2X,/S  (134) 

where  5  =  c{Tf-T^IL.  Note  that  the  quantity  e  =  k/(21iV^)  being  a  function  of  t  rather  than  q,  is 
the  source  of  nonsimilarity.  As  A  -» <»,  e  -»  0,  and  7(0,  t)  -»  T^,  and  the  problem  is  reduced  to  the 
classical  Stefan  solution  discussed  in  section  3. 1 ,  also  called  the  single-phase  Neumann  solution.  To 
solve  the  nonsimilar  eq  133,  134,  two-term  perturbation  expansions  are  assumed  for  F  and  X.  as 


F  =  Fo  =  eFi+0(£2)  (135) 

X  =  Xo  =  eXi +0(£2)  (136) 

Substituting  eq  1 35  and  1 36  into  1 33  and  1 34,  removing  the  implicitness  due  to  X  by  a  Taylor  series 
expansion  about  X^,  and  equating  the  coefficients  of  e°  and  e,  we  get 

Fo'+2qFo  =  0  (137) 

q  =  0,  Fo=r,  q=Xo,  Fo  =  0  (138) 

F,+  2qF|=0  (139) 
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Ti  =  0,  Fi  =  Fo;  Fi=- A,iFo(Xo) 

The  solution  of  eq  137  and  138  is  well  known  and  can  be  written  as 


(140) 


Fo=l-(erfTi/erfA.o)  (141) 

where  is  given  by  the  transcendental  equation 

Xocrf  X,oexp(x^)  =  S/V)i  (142) 

The  solution  for  F]  can  be  obtained  as 

F,  = - 2 - +  -l]  (143) 

Vn^erfXoL  erfXo  J 

where  X]  is  given  by  the  following  equation: 

Xi  = - ^ - .  (144) 

2e^Jc(erfXo)^+45(XoVjc^erfXo-tf^) 

In  order  to  compare  the  present  solution  with  the  solutions  shown  in  Figure  4,  we  calculated  X/e = 
/iXf/k  as  a  function  of  l/4£^  =  h^t/(kpc)  for  5= 0.1, 0.2, 0.5, 1 .0, 2.0  and  3.0.  TTiese  values,  shown  in 
Figure  4,  compare  well  with  other  approximate  solutions. 

4  J.  Freezing  of  a  semi-infinite  strip  (fin)  of 
liquid  not  initially  at  freezing  temperature 

Figure  6  shows  a  semi-infinite  strip  of  liquid  in  the  form  of  a  fin  with  cross-sectional  area  A  and 
perimeter  F.  The  liquid  is  initially  at  temperature  Tj.  At  time  r  ^  0  the  base  at  x  =  0  is  brought  and 
maintained  at  a  constant  temperature  that  is  lower  than  the  freezing  temperature  T(.  Both  the  solid 


Figure  6.  Geometry  for  the  freezing  of  a  liquid  fin. 


and  the  liquid  phases  convect  heat  to  the  surroundings  held  at  temperature  T^,  the  heat  transfer 
coefficient  being  h.  If  we  allow  for  the  surface  convective  heat  transfer,  the  transient  heat  conduction 
equations  for  the  solid  and  liquid  phases  can  be  written  in  partially  dimensionless  form  as 

dx^  ksA  tts  3r 


^ k0,)  =  A.<  x<  oo,  f  >0  (146) 

3x2  r  3f 
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where  Os  =  (7-s-  T{)/(Tf-  To),  0/  =*/  (T/-  Tf)/mT(-  To)),  Ba  =  {T»-  Tf)/{T{-  To) 
The  initial  and  boundary  conditions  are 


;rf(0)  =  0.  QiU0)  =  Bi=kt  {Ti-Tf)j{kATf-To)) 

6,(0, t)  =  -i,e,(x,t)= 0.  61  (X ,t) = 0 


dx  dx  vf-7’(,)is  dt 


Introducing  the  similarity  variables  t],  and  ti/  as 


(147) 


Tis  =  xli4a,t ,  T)/  =xTl4att  (148) 

and  defining  a  dimensionless  coordinate  (time)  e  =  4h  Pa,  t/(k,A),  we  seek  two-term  perturbation 
solutions  for  6$,  6/,  and  A.  as 

es(j:,r)=  i<o(ils)+£ui(t1s)  +  o(e^)  (149) 

0rW=  Vo(tl/)  +  evi(Ti/)  +  0(£2)  (150) 

X.(jc.r)=  V4as/  [Xo  +  eX,,+o(e2)]  (151) 

Substituting  eq  149  into  145  and  eq  150  into  146  and  removing  the  implicitness  in  the  interface 
boundary  conditions  by  Taylor  expansion  about  Xj,,  we  obtain  the  following  equations  for  Vp 

«o'+2t1s«;  =  0  (152) 

«o(o)=-l-  «o(Xo)  =  0  (153) 

'’o'+ 2ti^  v;  =  0  (154) 

Vo(aXo)=0,  Vo(«»)  =  9i  (155) 

UoiK) -  avp (a Xo)  =  2 pXo  (156) 

where  a  =  (asAx/)'^  and  P  =  L/(Cs(7'f-rp)).  Similarly  the  equations  for«|  and  V|  take  the  foim 
M|'+  2qs«,'-4«i  =  Uo-6a  (157) 

M,(0)  =  0,  M,(Xo)+ X,«p'(Xo)  =  0  (158) 

v,'+  2tl/  V,'-  4v,  =  Y  Vq-  6Ja^  (159) 

V I  (oo)  is  finite,  V|  (aXj  +  aXi  (aXo)  =  0  (160) 


(161) 


^  1  «o1^o )  +  )  -  “^^1 V  0  (a^o)  =  6  px, 

where  y=q/Cs. 

Equations  152-156  constitute  the  classical  two-region  Neumann  problem  whose  solution  is 


erf 


(162) 


Vo  =  9i-ei-^^p^-  (163) 

erfc  ^01  ^ol 


where  Xq  is  given  by  the  transcendental  equation 


e”X>/erf  X,o-0i  a e““  X)/erfc  (aXo)  =  1/71^3X0  (164) 

Churchill  and  Evans  (1971)  give  the  solution  of  eq  164  for  a  range  of  the  parameters  involved. 

The  solutions  for  Mj  and  V|  appear  in  the  form  of  repeated  integrals  of  the  error  function  and  are 


Ml  =fli/2(qs)-(l +0a)  /^erfe^s-i  -(erfcqs)/4erfXc 

L  4, 


=  Ci/^erfcTi^  +9a/(4a2)_y0j[i  _erfqif/erfc{aXo)]/4 


(165) 

(166) 


where 

In  (tIs)  =  '■ "  erfc  (-Tls)  +  (- 1 ‘  "  erfc  . 


=  [( l+0a)<^  erfc Xo- 2X.oX,  -0a/4]//2(X.„), 

C,  =  [-2(>lo-PXo)  X,  -0a/(4a^)]//2erfc(aXo) 

The  solution  for  X]  takes  the  form 

Xl  =  {[(l+  0a)^o]/[l+  2Xo(Xo+  -4o)] 

-0a[7,  +A-,/a2]/4-t-Y{>l„-pxJ/2-/lo/2}/L,  (167) 

where 

/\o  =  e“W(Vn  erf  Xo),  J\  =1\  (Xo)//2(Xo) 

Ki  =  a/erfc(aXj//^erfc(aXo) 

L,=6p  +  4(l-a2)Xo/lo  +  X4oUi+^i)  +  2(2a2Xo-X:,)pXo 
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It  is  interesting  to  note  that  unlike  the  Neumann  problem,  the  freezing  front  can  become  stationary 
if  >  Tf.  This  would  occur  when  the  convective  heat  transport  from  the  surroundings  to  the  fin  from 
the  top  and  bottom  faces  matches  the  heat  extracted  through  the  base  of  the  fin.  Under  the  condition 
of  a  stationary  freezing  front,  the  temperature  distributions  through  the  solid  and  the  liquid  phases  are 
given  by 

0s=0a-{l  +  0a)coshpj:+[(l+ea)coshfiA.„-0a](sinhpj:)/(sinh|iA,)  (168) 

0/  =  [l  -  (169) 

where  =  (hP/kgAY^,  K  =  kjki  and  (the  stationary  location  of  the  freezing  front)  is  given  by 
the  explicit  relation 

A,^={l/fi)ln[c+(c2_),)i'2]  (170) 

where  b  =  (Ar-1)/(A:+1)  and  c  =  A:(l+0a)/((l+AO6a). 

Figure  7  shows  a  typical  result  for  the  progress  of  the  freezing  front  for  a  =  0.7 1 ,  P  =  2.5  (Stefan 
number = 0.4),  /( = 0.7 1 , 0j  =  1 ,  and  % = 0.5.  This  case  corresponds  to  a  solder  fin  for  which  Tj, = 82®C, 
Tf  =  183°C,  and  T;  =  Tf^  =  233.5°C.  Examination  of  Figure  7  reveals  that  the  two-term  perturbation 
solution  virtually  coincides  with  the  finite  difference  solution  up  to  e  =  4.  However,  as  e  increases 
further,  the  perturbation  solution  begins  to  deviate  significantly  from  the  finite  difference  solution.  At 
e  =  10,  the  error  is  about  13%.  The  zero-order  solution  (Neumann)  is  also  shown  for  comparison  and 
it  is  evident  that  the  first-order  correction  improves  the  accuracy  considerably.  Figure  8  shows  the 
corresponding  temperature  profiles.  Again  the  perturbation  solution  is  quite  accurate,  up  to  e  =  3,  but 
the  accuracy  deteriorates  considerably  at  e  =  1 0,  particularly  in  the  liquid  region.  Figures  7  and  8  have 
been  adapted  from  Rathjen  and  Jiji  (1970),  who  also  provide  a  graph  of  base  heat  flux  vs.  e. 


Figure  7.  Comparison  of  perturbation  and  numerical  solu¬ 
tion  for  fin  freezing. 
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Figures.  Temperature  profiles  during 
freezing  of  fins. 


5.  SINGULAR  PERTURBATION  EXPANSIONS 


The  term  singular  perturbation  expansion  is  used  to  describe  expansions  that  lack  the  feature  of 
uniform  validity  thatcharacterizes  the  regular perturbationexpansion.  In  sections  3  and 4,  the  solutions 
obtained  were  valid  throughout  the  domain  of  the  independent  variable.  Now  we  study  problems  for 
which  the  expansion  fails  in  certain  regions,  called  “regions  of  nonuniformity"  m-  “boundary  layers.” 

The  nonuniformity  exhibits  itself  in  several  forms:  1)  the  solution  becomes  infinite  at  some  value 
of  the  independent  variable,  2)  the  solution  has  a  discontinuity  within  the  domain  of  interest,  3)  the 
sol  ution  fails  to  satisfy  some  boundary  condition,  and  4)  the  solution  contains  an  essential  singularity. 
Two  examples  of  nonuniformities  arising  in  freezing  problems  are  discussed  below.  Other  examples 
will  be  discussed  when  the  various  techniques  of  handling  singular  expansions  are  presented. 


5.1.  SoUdificatkMi  of  a  saturated  liquid  in  a  spherical  dcnnain 

Consider  a  saturated  liquid  of  inrinite  extent  outside  a  sphere  of  radius  as  shown  in  Figure  9. 
At  r  >  0  the  surface  at  is  suddenly  lowered  to  a  subfieezing  temperature  T^.  Fot  the  solid  phase, 
we  have 

=  (171) 

R  dR^  “ 
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T(«..l)=  T„,  r(«,.0»  T,.  k  |I|  =  (>t 

Introducing  the  following  dimensionless  quantities 

(172) 

«={r-ro)/(rf-ro),  r=«//?^.  x  =  k{Tf-T^)t/pLRl 

into  eq  17 1  and  172  and  replacing  the  variable  r  by  the  variable  Rf,  we  get 

(173) 

I  d\ur)  _  j.  3m  du 
f  dr^  drt  dr 

(174) 

«(r=l,/f)  =  0,  M(r=rf,rf)=l. 

Assuming  a  two-term  perturbation  expansion  of  the  form 

(175) 

«=  Mo+  £mi  +  £^“2 

(176) 

and  substituting  into  eq  174,  175,  the  governing  equations  for  ug,  U]  and  U2, 

can  be  obtained  as 

o 

ti 

»— f|  k 

(177) 

Uo{r=  1 ,  rf)  =  0 ,  «o(f=  /t)  =1 

(178) 

1  _  3«o  3«o 

dr^  dri  dr 

'  r-r\ 

(179) 

«!('•=  1 . '•f)  =  0,  M,(r=tj, /t)=0 

(180) 

1  3^(«2r)  3io  .  3«i  3«o 

r  dP-  ~  dvi  '  dr  ■  3r 

'  r=rf  ‘  r=r\ 

(181) 

U2{>'=  1 ,  rir)  =  0,  U2(r=  rt,  r[)  =  0 

(182) 

Solving  eq  1 77-1 82,  the  solution  for  u  is  obtained  as 


(183) 


Using  eq  1 83  to  evaluate  and  then  integrating  eq  175,  the  solution  for  x  can  be  obtained  as 
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T  =  ^'^.  +1  (184) 

6  6  45  »r 

where  the  initial  condition  rf  =  I ,  x  =  0  has  been  utilized. 

The  solutions  represented  by  eq  1 83  and  1 84  were  first  reported  by  Pedroso  and  Domoto  ( 1 973c). 
They  are  uniformly  valid  for  outward  spherical  solidification.  However,  for  inward  spherical 
solidification,  the  second  term  in  eq  1 83  becomes  singular  as  the  freezing  front  approaches  the  center, 
that  is,  as  0.  The  divergence  of  the  sol  ution  as  »  0  is  shown  in  Figure  1 0.  This  singularity  grows 

further  in  the  third  term  because  of  the  presence  of  the  term  rf  in  the  denominator.  Examination  of  eq 
1 84,  however,  shows  that  the  freezing  time  solution  is  uniformly  valid  up  to  (He)  and  the  singularity 
first  appears  in  the  third  term.  Both  eq  1 83  and  1 84  must  be  rendered  uniformly  valid  if  they  are  to  be 
used  for  complete  inward  freezing.  This  will  be  achieved  later  using  the  method  of  strained 
coordinates. 


Figure  10.  Freezing  radius  vs  time  for  inward  spherical 
freezing. 


52.  Solidification  of  a  saturated  liquid  in  cylindrical  domain 

The  problem  of  section  5. 1  is  now  considered  for  cylindrical  geometry.  The  temperature  distribu¬ 
tion  in  the  solid  phase  is  governed  by 

=  (185) 

R  dR\  dR  I  a  a, 

with  boundary  conditions  dictated  by  eq  1 72.  Using  die  dimensionless  quantities  eq  1 73,  the  problem 
is  reduced  to 


lA 

r  dr 


(ry 

=  e—  . 

du 

1  dri 

d/r 

dr 

(186) 


with  boundary  conditions  established  by  eq  1 75.  Repeating  the  procedure  of  the  previous  section,  the 
two-term  perturbation  solution  for  u  is  obtained  as 
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„  =  inL +  e 


Inrf 


1+  r^(ln/ir-l ) 

4r^{lnrf)'‘ 

Ar^(lnrf)^  I 

(187) 


If  we  use  eq  1 87  to  calculate  ^ 


and  then  integrate  eq  1 75,  the  solution  for  T  subject  to  the  ini- 


r=rt 


tial  condition,  rf  =  1,  x  =  0,  is  obtained  as 


x  =  l-r^  In/y  +  J-  (l-  rf^)+  —  e(l-  r^)(l  +  — L-\  (188) 

2  4  4  \  Inff/ 

Equations  187  and  188  are  unifoimly  valid  for  outward  cylindrical  solidification.  These  same 
equations  are  valid  for  the  inward  cylindrical  solidification  but  they  break  down  near  the  end  of  the 
solidification  process  because  if  -» 0,  the  second  term  in  eq  1 87  becomes  singular.  If  one  calculates 
the  third  term  of  the  expansion,  singular  behavior  as  »  0  worsens  because  of  the  appearance  of  the 
higher  powers  of  rf  in  the  denominator.  The  situation  parallels  the  one  discussed  in  the  previous  section 
for  the  inward  spherical  solidification.  The  two-term  solution  for  x  given  by  eq  1 88,  however,  retains 
its  validity  for  complete  inward  cylindrical  solidification,  although  if  the  third  term  is  added  to  eq  1 88, 
it  exhibits  singular  behavior  near  the  end  of  the  solidification  process.  We  will  return  to  the  present 
problem  in  the  section  on  the  method  of  strained  coordinates  where  uniformly  valid  solutions  are 
developed  for  the  inward  case. 


6.  MISCELLANEOUS  EXAMPLES  OF 
PERTURBATION  EXPANSIONS 

Besides  the  problems  discussed  so  far,  perturbation  methods  have  also  been  used  for  other  phase 
change  problems.  This  section  is  a  review  of  such  studies  with  details  in  many  cases  left  to  the 
appropriate  references.  However,  some  of  the  problems  will  be  examined  in  detail  in  the  sections  that 
follow. 

For  one-dimensional  freezing  in  planar  geometry,  Huang  and  Shih  (1975)  developed  a  three-term 
perturbation  solution  for  the  freezing  of  a  warm  fluid  over  a  flat  plate  cooled  from  below,  and  found 
that  the  growth  of  the  solid  phase  predicted  by  the  perturbation  analysis  agrees  with  the  experimental 
results  of  Siegel  and  Savino  (1966)  for  time  r  >  1(X)  seconds,  but  for  f  <  1 00  seconds  the  perturbation 
solution  overestimates  the  growth.  A  variation  of  the  problem  discussed  in  section  3.5  is  the  subject 
of  a  paper  by  Gutman  (1986).  He  considered  the  freezing  of  a  finite  slab  cooled  on  both  sides  by  a 
coolant  whose  temperature  changes  sinusoidally  with  time.  Because  he  assumed  that  the  freezing  is 
independent  on  each  side  of  the  slab,  he  was  able  to  demonstrate  that  his  solution  for  the  case  when 
the  far  end  of  the  slab  is  insulated  matches  with  the  two-term  expansion  of  the  exact  Stefan  solution 
(section  3.1).  Gutman  also  considered  the  case  of  a  finite  slab  with  constant  heat  flux  at  one  face  while 
the  other  face  is  insulated.  The  perturbation  approach  used  in  section  3.3  can  also  be  used  for  other 
surface  temperature-time  variations.  For  example.  Lock  (1969)  gives  a  three-term  solution  for  the 
power  law  variation,  T(0,r)  =  Cr",  where  C  and  n  are  constants.  Indeed,  he  showed  in  another  paper 
that  his  method  holds  for  arbitrary  variations  of  the  surface  temperature  (Lock  1971). 

The  analysis  for  spherical  and  cylindrical  solidification  with  constant  surface  temperature, 
presented  in  sections  5. 1  and  5.2,  has  been  extended  to  the  case  of  convective  cooling  of  the  wall  by 
Huang  and  Shih  (1975).  They  found  that  for  outward  growth  of  the  solid  phase,  the  solutions  are 
uniformly  valid  but  the  same  solutions  when  applied  to  the  inward  freezing  situation  begin  to  diverge 
towards  the  end  of  the  freezing  process.  The  case  of  a  more  general  boundary  condition  involving 
simultaneous  convective  and  radiative  cooling  (or  heating)  has  been  considered  by  Seeniraj  and  Bose 
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(1982).  However,  the  sol  utions  for  both  cylindrical  and  spherical  geometries  are  restricted  to  only  two 
terms. 

The  difficulty  associated  with  the  singular  behavior  of  the  regular  perturbation  solutions  for  inward 
phase  change  problems  has  been  addressed  by  several  authors.  For  example,  the  singularity  of  eq  1 83 
and  1 84  for  inward  spherical  freezing  has  been  remedied  by  Pedroso  and  Domoto  (1973d)  with  the  aid 
of  the  method  of  strained  coordinates.  The  same  method  was  used  by  Milanez  and  Boldrini  ( 1988)  to 
obtain  a  uniformly  valid  solution  for  inward  freezing  of  a  sphere  with  convective  wall  cooling.  For 
inward  cylindrical  freezing  with  constant  temperature,  the  singular  solution  given  by  eq  187  was 
rendered  uniformly  valid  by  Afsar  et  al.  (1979)  who  developed  an  appropriate  strained  coordinates 
solution.  In  a  more  recent  paper.  Parang  et  al.  (1990)  have  reported  strained  coordinate  solutions  for 
both  inward  cylindrical  and  spherical  solidification  when  the  freezing  at  the  wall  is  accomplished  by 
simultaneous  convection  and  radiation.  A  unified  approach  for  inward  solidification  that  allows 
simultaneous  treatment  of  the  problem  in  plane,  cylindrical,  and  spherical  geometries  with  boundary 
conditions  of  constant  temperature,  constant  heat  flux,  or  pure  convection  was  developed  by 
Prud’homme  et  al.  (1989).  Again,  the  tool  used  to  derive  uniformly  valid  solutions  was  the  method  of 
strained  coordinates. 

Another  technique  that  has  been  employed  to  treat  singular  perturbation  expansions  is  the  method 
of  matched  asymptotic  expansions.  Weinbaum  and  Jiji  (1977)  used  this  technique  to  study  the  freezing 
of  a  finite  slab  not  initially  at  the  freezing  temperature.  The  boundary  conditions  chosen  were  those 
of  a  constant  subfireezing  temperature  at  one  face  while  the  other  face  was  either  insulated  or  kept  at 
the  initial  temperature.  Since  even  the  first-order  term  of  Weinbaum  and  Jiji’s  solution  was 
complicated,  Charach  and  Zoglin  (1985)  approached  the  same  problem  by  first  developing  the  heat 
balance  integral  formulation  and  then  constructing  a  perturbation  expansion.  With  this  strategy  they 
were  able  to  determine  higher-order  terms  more  easily.  Gutman  (1987)  attacked  the  same  problem 
using  the  classical  Neumann  solution  for  the  two-region  problem  as  the  inner  solution,  and  a  modified 
solution  as  the  outer  solution.  The  modification  to  the  Stefan  solution  involved  the  replacement  of  the 
time  variable  r  by  (r  -  C)  where  the  constant  C  accounted  for  the  additional  time  needed  for  complete 
freezing  because  of  initial  superheating  of  the  liquid,  that  is,  Tj  >  Tf .  By  matching  the  inner  and  the 
outer  solutions  using  the  overall  energy  balance  as  a  criterion,  Gutman  was  able  to  find  the  constant 
C.  In  another  study,  Jiji  and  Weinbaum  ( 1 978)also  used  the  method  of  matched  asymptotic  expansions 
for  freezing  in  an  annular  region  with  liquid  not  initially  at  the  freezing  temperature.  A  constant 
temperature  was  imposed  at  the  outer  surface  while  the  inner  surface  was  assumed  to  be  either 
isothermal  or  adiabatic.  However,  unlike  their  earlier  analysis  for  the  plane  geometry,  the  method 
failed  to  give  a  uniformly  valid  solution  for  the  annular  geometry.  The  solution  behaved  well  initially 
but  as  the  freezing  front  approached  the  inner  surface,  that  is,  as  the  last  bit  of  liquid  froze,  a  singularity 
appeared.  Such  behavior  is  to  be  anticipated  because  the  liquid  phase  degenerates  at  this  point. 

The  method  of  matched  asymptotic  expansions  has  also  proved  effective  in  dealing  with  the 
singular  nature  of  the  perturbation  solution  for  inward  spherical  solidification.  For  example,  Riley  et 
al.  ( 1 974)  corrected  eq  1 83  and  1 84  in  a  multi-region  structure  using  this  method.  However,  Stewartson 
and  Waechter  (1976)  later  found  that  the  inner  expansion  developed  by  Riley  et  al  :.lso  fails  in  a  minute 
region  just  before  the  center  freezes  and  must  be  supplemented  by  an  additional  expansion  in  that 
region.  Thus  a  completely  valid  solution  consists  of  three  expansions  rather  than  two.  4  similar 
approach,  but  unified  for  both  spheres  and  cylinders,  has  also  been  presented  by  Soward  (1980). 
Howarth  (1987)  adopted  the  Riley  et  al .  ( 1 974)  method  to  study  inward  spherical  freezing  under  the 
condition  of  constant  heat  flux  and  reached  the  conclusion  that  the  breakdown  of  the  outer  solution  for 
the  constant  heat  flux  condition  occurs  at  0  (e*^'*)  rather  than  0  (e*^)  as  is  the  case  for  the  constant 
temperature  condition. 

Another  class  of  phase  change  problems  where  the  perturbation  techniques  have  proved  effective 
is  the  two-dimensional  Stefan  problems.  The  basic  idea  is  to  develop  a  two-dimensional  solution  as 
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a  perturbation  of  a  one-dimensional  case.  For  example,  Schulze  et  al .( 1 983)  have  shown  how  the  two- 
region  planar  freezing  of  a  liquid  with  slightly  varying  wall  temperature  or  heat  flux  can  be  treated  as 
a  perturbation  of  the  Neumann  solution.  Similarly,  Howarth  (1990)  considered  freezing  on  a  slightly 
wavy  wall  as  a  linearized  perturbation  of  the  Neumann  solution.  The  problem  of  two-dimensional 
solidification  in  a  comer  can  also  be  analyzed  as  a  perturbation  of  the  Neumann  solution  (Howarth 
1985). 

When  the  liquid  flows  over  the  cold  surface,  the  phase  change  process  is  controlled  simultaneously 
by  conduction  and  convection.  Despite  the  complexity.  Lock  and  Nyren  (1971)  found  that  a  regular 
perturbation  approach  is  still  useful.  In  particular,  they  considered  the  fully  developed  flow  of  a 
freezing  liquid  in  a  long  circular  tube  cooled  by  external  convection.  The  perturbation  parameter 
was  chosen  as  e  =  SBil{  1  +Bi),  where  S  is  the  Stefan  number  and  Bi  is  the  Biot  number  representing 
the  external  convection.  Closed  form  solutions  are  given  for  the  zero-order  and  the  first-order 
problems.  The  growth  of  solid  phase  due  to  sudden  lowering  of  the  wall  temperature  in  fully  developed 
laminar  flow  of  a  liquid  in  a  circular  tube  has  been  studied  by  Cervantes  et  al.  (1990)  using 
kJkt{Tf-Tj{To-T^)  as  a  perturbation  parameter,  where  represents  the  initial  temperature  of 
the  liquid  and  the  other  symbols  have  their  usual  meaning.  They  generated  a  two-term  solution  which 
is  valid  when  the  thickness  of  the  frozen  layer  is  small  compared  to  the  tube  radius. 

To  conclude  this  section,  we  finally  refer  to  the  applications  of  perturbation  methods  to  soldering 
and  welding  problems.  Andrews  and  Atthey  (1975)  model  the  process  of  hole  formation  by  a  laser  or 
an  electron  beam  as  a  planar  evaporating  boundary  heated  by  a  constant  power  density  source.  The 
model  was  solved  as  a  two-term  perturbation  expansion  in  which  the  ratio  of  sensible  heat  to  raise  the 
material  to  the  evaporation  temperature  and  the  latent  heat  is  taken  as  a  perturbation  parameter.  The 
method  of  matched  asymptotic  expansion  was  used  to  obtain  a  uniformly  valid  solution  for  the  velocity 
of  the  evaporation  front.  A  similar  technique  has  been  used  by  Antaki  (1990)  to  predict  the  laser- 
induced  sublimation  of  a  finite  layer  of  solid  material. 


7.  METHOD  OF  STRAINED  COORDINATES 

The  basic  idea  underlying  the  technique  is  to  expand  both  the  independent  and  dependent  variables 
in  terms  of  e  with  coefficients  expressed  as  functions  of  a  new  independent  variable.  The  coefficients 
of  the  independent  variable  series  are  called  the  straining  functions.  The  expansions  are  next 
substituted  into  the  original  equations  to  generate  the  usual  sequence  of  perturbation  equations.  It  is 
at  this  stage  that  the  choice  of  the  straining  functions  is  made  such  that  higher  approximations  are  no 
more  singular  than  the  first.  This  principle  is  often  referred  to  as  Lighthill’s  rule,  and  its  application 
is  the  crucial  step  of  the  whole  analysis.  If  successful,  the  result  is  an  implicit  but  uniformly  valid 
solution.  Because  of  its  first  appearance  in  Lighthill’s  paper  (1949),  the  method  is  also  called 
Lighthill’s  technique. 

The  spirit  of  Lighthill’s  technique  is  also  reflected  in  earlier  works  of  Lindstedt  ( 1 882)  and  Poincard 
( 1 892)  where,  instead  of  a  coordinate,  a  parameter  is  strained  to  achieve  uniform  validity.  When  a 
parameter  is  strained,  the  technique  may  be  appropriately  termed  the  method  of  strained  parameters. 
Giving  credit  to  the  contributions  of  Poincar^ ,  Lighthil),  and  later  works  of  Kuo  ( 1 953, 1 956),  the  name 
PLK  method  was  coined  by  Tsien  (1956). 

We  illustrate  the  application  of  the  technique  to  several  problems.  The  first  problem  is  discussed 
in  full  detail  but,  for  the  remaining  ones,  the  essential  steps  are  indicated  with  discussion  of  results. 

7,1.  Inward  spherical  solidification  with 
constant  wall  temperature 

In  section  5. 1 ,  a  regular  perturbation  solution  was  developed  for  outward  spherical  solidification, 
and  it  was  noted  that  the  .same  solution  applies  to  inward  solidification,  but  the  solution  becomes 
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singular  as  the  center  freezes.  This  deficiency  can  be  remedied  by  straining  the  coordinates  r  and  ff, 
together  with  the  dependent  variable  u.  Following  Pedroso  and  Dcnnoto  (1973d),  we  introduce  two 
new  variables  (|)  and  and  expand  u,  r,  and  rf  as  follows: 


«=  r/o(^,y)+  ewi  (<|),y)+  e2M2(<)),y) 

(189) 

r=<|>+  eCTi(<|>,y)+  e2a2{<|i,y) 

(190) 

rf=y  +  eai(y,y)+  £2o2(y,y) 

(191) 

where  the  straining  functions  ai(<t),y)  and  C2(<t>,\(r)  will  be  chosen  to  ensure  uniform  validity  of  the 
solution.  Note  that  as  r  — » /-f,  <|>  — >  t|r.  Since  the  regular  perturbation  solution  is  satisfactory  near  the 
spherical  boundary,  we  choose  to  make  the  straining  vanish  at  the  boundary  so  that 


lim  ai{4>,t|/)  =  0  (r-»l,<|)->l) 


(192) 


lim  ai(\|/.v)  =  0  (rf^l.\)/^l) 

Also  note  that  r  is  a  function  of  both  <|)  and  y  but  rf  is  a  function  of  y  alone. 

To  change  the  variables  from  (r,  rf)  to  (<t>,  y)  in  eq  174  and  175  we  proceed  as  follows.  For  any 
function /(r,  rf)  where  r  =/i(<t>,  y)  and  rf  =  /2(y),  we  have 


3/  _  ^  ^  ^ 

3(t>  3r  3<t>  3rf  3(t) 


(194) 


3y  3r  dy  3rf  3y 
Since  rf  =/2(y),  rf/3<t>  =  0;  hence 


(195) 


3/  _  ^  ^ 

3(t)  3r  3<t) 

From  eq  196  and  195  we  have 


3/_  3/  1 3r\  * 
3r  3(t)  \3(t)/ 


or 


3/t 


3rf 


3/  3/  3r|  /3rfV 

3y  3r  3y/  v3y/ 


1-1 


3f  3r 


3y  3<t)  \3(|>/  3y 


-'_^1  /^r' 

3yJ  \3y/ 


(196) 


(197) 


(198) 


Since  we  need  the  second  derivative  with  respect  to  r.  we  differentiate  eq  197  to  get 
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ar2 


ay  9/a^rj 

liL] 

a<t>^  a<i>  a(|)^ ' 

la(t>/ 

(— r 


(199) 


Foreq  174  the  quantities  required  are  du/dr,  du/drf,  and  d(ur)^r^.  Based  on  eq  197-199,  we  can  write 


du . 
dr 


(200'* 


du  ra«  dufdr\-^dr 


d/t  L5V 


(i 


(201) 


_ 

dp- 


di^iru)  d(ru) 


[  3<t>  '3v 


(202) 


Using  eq  189-191,  we  derive  the  following  derivatives 


du  9 

- - M0^+  £'‘U2^ 

a<ii 

(203) 

U0ur+  e«iw+ 

av 

(204) 

— =  1  + 

<79 

(205) 

ea,^  +  e2a2,t,,|, 

d<p 

(206) 

— =  eoiv  +  £^2\u 

av 

(207) 

1  +  eoiy  +  e2a2y 

(208) 

where  the  subscripts  <ti  and  are  used  from  now  on  to  indicate  partial  derivatives,  and 

a,  =  ai(<ti,v) 

02=02(4*,^) 

(209) 

CTi  =  a,(t|/,v) 

02  =  a2(v.V) 

Additionally,  we  need  the  derivatives  d(ru)ld^  and  d\ru)ldfy^.  From  eq  189  and  190,  we  have 
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ni=  (^  +  eoi  +  e2o2)(i4o+  e«i  +  t^u2) 


or 


rus(^iiQ+  8(uoOi  +  ♦«!)+  e^(Mo02+  “1®1+  ♦“2) 


Thus 


9(>'«)- 

d<|> 


=  (<|>Uo)^+  e(uoOi  +  <|>Ui)^+  £2(11002+  MiOi  +  (>U2)^ 


(210) 


and 


3^(ru)_ 
a<t»2 

Returning  to  eq  200,  we  have 


=  (<t>«o)^+  e(«o<Ji  +  e2(Mo02+  uiOi  +  (|>M2)^ 


(211) 


or 


^={"0*+  £“!♦+  £^«2+)(1  +  £«!♦+  ‘ 

dr 


du 


au  /  \ 

— =Uo^+e(H,^-«o^a,^) 


(212) 


where  we  have  retained  only  ternis  up  to  (Ke).  Sincee  appears  on  the  right-hand  side  of  eq  1 74,  we  need 
to  retain  only  two  terms  in  eq  212  to  obtain  the  solution  up  to  0(e2). 

We  now  work  out  the  expansion  for  du/dr {.  Using  eq  201 ,  and  retaining  terms  only  up  to  0(E),  we 
get 


or 


du 

drt 


du  (  ^ 

**oy'*’  ^V“iy~  “ov®iV' 

drt 

Finally,  coming  to  eq  202,  we  have 

{('l'«o)^+  e(«oO!  +  <t>«i)^+  e2(Ho02+  uioi  +  ^uz)^ 

-[(<t‘“o)^+  e{«0Oi  +  4n<i)^+  e^(Ho02+  “lOi  +  ^«2)J 
x(eo,^+  e2a2^)[  1  -  eai^+  02^)] ) 

X  [  1  -  e2oi^  +  £2  (30|^  -  2<T2^)] 


=  [(«oy  +  e“iy)-  («o4,+  e«i^) (1  -  eo^)Eaiy]  (l  -  £Oiy) 


(213) 
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or 


Oi«^{<l'“o)^-  2oi^(<|>«oy 


+  e2[^(<t>«2  +  MiOi  +  «0<J2),(i4,-  <y|#(<l>“l  +  “OOl)^ 

-2oi(|.(<t>«i  +  «oOi)^-  (2024,-3  of^)((t>uo)^ 

-  (024*1,  “  <ri4iiti)(<t>«o)J 

The  right-hand  side  of  eq  174,  when  multiplied  by  r,  can  now  be  written  as 


(213a) 


£®l)[«04>l4»=v+  e(«14>-  “04.Ol4.)l4.=y] 

v/"  r=rf 

X  [l/Qx^  +  6(ui\),—  no^O\y—  MoyOjy)] 

=  e  {<l>Wo^l4,=  \j,UO\(f+  £[^^)(>I$=4»(*^1\|»~ 

+  Mov(®l“0<>l4)=v+  (“14~  “04iOl4>)  '4i=v*t*)]} 

We  now  use  eq  213a  and  214  to  write  the  sequence  of  perturbation  equations  as  follows: 
e°:  ('l'«o)^=0 


du 


(214) 


(215) 

(216) 

(217) 

(218) 


«o(v.4'=l)=0  Mo(v.<t>=  v)=l 

e' :  {^\  +  mqOi)^-  Oi44,(<1)«o)4,=  <|>M04,l4*=vMoy 

Ml(v,(|>=l)=0  Mi(\|r,<|»=  v)  =  0 

e^:  (^2  +  «iO]  +  «o02y-  Oi44,(<)>«i  “oOiy  2oi4,(<|>mi  + 

-(02#-  30l,jOl,^)((t>Mo)4=  <t>«04l^=y(Mly-  M0^O]y 

-«Ov®lv)+  «0v[oi«04)l4i=v+  <t>(“l<>“«0^<7l4)4,=y]  (219) 

“2(¥><t>=  l)=(J  M2(¥.<I>=  ¥)  =  0  (220) 

Zero-order  solution.  Since  eq  215  and  216  have  the  same  form  as  eq  177  and  178,  the  solution  is 

„0=i:^  (221) 

l-l/\|t 
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First-order  solution.  From  eq  221  we  have 


«o^  = 


_ }1L_ 


1 

V(l-  ¥) 


<t>M0  = 


¥(<l>-l) 

\|f-l 


<H  =- 


* 


!-<!> 


Using  eq  221  and  222,  eq  217  can  be  simplified  to 


+_!t_5L\  =_ 

1-  ¥  <l>/^ 


!-<!> 

¥(l-¥f 


Integrating  eq  223,  we  have 


- 

\  1-  y(l-  yf 


Let  u\  be  identically  zero.  Then,  imposing  the  condition 


A(¥)  = 


1 


Integrating  once  again,  we  have 


^*™j(cti/<l>)4,  =  0,  gives 


- L+  /2(y) 

1  -  ¥  <1>  6\|/(  1  -  2\|>(  1  - 

Using  the  condition  in  eq  192  gives 


^2(¥)  = - i - - 

6¥(  1  -  ¥P 

Hence  the  solution  for  Oj  becomes 

c.-JtOjJlL 

Second-order  solution.  By  making  U2  identically  zero,  eq  219  reduces  to 

Oi#(«oOi)^-  2oi,|,(Moai)^- {o20(j,-  3ai^oi^)((l)Mo)^ 

=  -  l|>U0^l,^=y(M0^Oiy+  M0yO|y)+  HOy[oi  “04>l^=y 


(222) 


(223) 


(224) 


(225) 
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We  now  evaluate  the  various  quantities  appearing  in  eq  225.  From  eq  224,  we  have 


^  _(l-»)"(4i|>-l). 


(l-»)(l-2») 

y2(i_y)2 


From  eq  221  and  224  it  follows  that 


«0Ol=- 


6\|/(l-v)^ 


and 


- 

*  3v(l-v)^ 


:  (“oOi)^ 


V(l-¥P 


Also,  from  eq  224  we  obtain 


Oi  =  ai(¥¥)=-^^-^ 

6v 

^  _<t>(l-<t>r(l-2x|f) 
- 

3v^(l-Vp 


and  Oiy  = 


6¥^ 


Using  the  foregoing  information  we  can  evaluate  the  various  terms  appearing  in  eq  225  as  follows 

*  3v’(l-v)’ 


2\|r^(l-v)^ 

.  .  /  -  \  {2-4w)(l-<bp  1-(|) 

<fr“o  i|)=v(«o<|i®iv''  - — - - — - — 

6\|rMl-¥r  6v^(l-vr 

«0v  (oi  --37^-73-  • 

6y^(l-¥r  6v^(l-vr 

Making  use  of  the  foregoing  information  and  simplifying,  eq  225  now  appears  as 
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V  (02\  _(4-.-4v-15»Xl-<l>P  ,  2(1- ♦) 

6v*(l-v)*  3v^(l-\|^ 

Integrating  eq  226,  we  get 

V  =  0+*f(l2*-4^^  _*(2z1L+ ft(« 

1  -  v  I  <t»  24\|»*(  1  -  3t|^{  1  -  vf 


Imposing  the  condition,  (02/4*)^  =  0,  gives 
^  1 


/3('l0  =  - 


3t|»2(2-  v)3  ■ 


Integrating  once  again,  we  get 

V  <32  -  (1-  <|>)^(l-4y+10(|>)  ^  _ <l>  |/4(W1. 

1-V4>  120\(^(  1  -  9v^{l-V)^  3v^(l-Vp 

Using  the  condition  eq  192  gives 


V{i- v)^ 

Thus  the  solution  for  C2  becomes 

o,  =  fi-  (1-«|>P(1-4V+10«|>) 

"\|P(1-V)2[9  120\|»(l-v)^ 


From  eq  227  it  follows  immediately  that 


02  =  a2(V.'rt  = 


(22v-3)(l-v) 

360\|i^ 


Also,  for  subsequent  use,  we  deduce 


(227) 


^  ,  -22a|/^  +  10w-3.  ^  22y^-50w+9 

<32*I*=v= - ^ - .  <32v= —  - T - • 

^  360^*  360V* 

Having  determined  the  straining  functions,  we  can  now  write  the  final,  formal  solution  from  eq 
189-191  as 


„=llli^  +  0(e3) 
1-1/V 


r=  e 


4>(1-4>P 


+ 


1  (1- ♦l^ll +4v+ 104») 

N^(i- vy' 

9  120v(l-t|^ 

(228) 


(229) 
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(230) 


I 


rf-  y-  F  -L-  e2  (22v-3  )( 1  -  y) 

6v  360V 


Freezing  time.  To  obtain  the  freezing  time  x  as  a  function  of  \|/,  we  consider  the  third  condition  in 
eq  175.  Since 


dz  ek^\(A\fj 
we  have 


dL=i^  ’  £!1 

{drr=rfl 


(23i) 


To  obtain  the  expansion  for(du/d/lr=rf)  up  to  0(£^),  we  return  to  eq  200.  Noting  that  and  M2  ^re 

identically  zero,  we  have 

^=uoip(l+  eOi(>+  e2o2,^)~* 
ar 

=  «o« -  eoi(t.M04. -  e^[ (<J2« -  o? Jko,,,]  +  0  (e3)  (232) 

Hence 

=  tto<,lo=v-  <^33) 

r=rf  ^  ^ 

Substituting  eq  233  into  eq  23 1  we  have 


du 

dr 


{«o.frV=v-  J*’  (l  +  ec,y+  e2S2y) 

which  can  be  expanded  and  simplified  to  give 


^=(“o^4=v)"‘[i+  £(<3l^l«=V+  Siv) 


+  e2  (o2,^+  <Ti^l(),=  y<Tn(r+  C2^l^=y)] 

Utilizing  the  appropriate  expressions  for  quantities  appearing  in  eq  234,  it  is  found  that 

^  =  -\|/(l-v)-e2-(i  -  y)  +  e2j_/j - L\ 

d\^  3  90\rjr2  ^J^2| 


Integrating  eq  235  and  imposing  the  condition  y  =  1,  x  =  0,  the  final  result  is 


6  3  180^2 


(234) 


(235) 


(236) 
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Figure  11.  Temperature  distribution 
for  inward  spherical  solidification. 


Figure  12.  Freezing  radius  vs  time  for 
inward  spherical  solidification. 


Equations  22S-230,  and  236  constitute  the  complete  uniformly  valid  solution.  Since  the  solution  is 
implicit,  the  computation  proceeds  as  follows.  First,  values  of  e  and  rf  are  fixed.  Next,  eq  230  is  solved 
by  iteration  to  obtain  the  value  of  \|t.  Choosing  values  of  <|>  in  the  range  t|/  to  1 ,  eq  229  and  228  are  used 
to  calculate  the  temperature  distribution  u  as  a  function  of  r.  The  freezing  time  is  calculated  using  eq  236. 

Comparison  with  numerical  solution.  Sample  results  for  the  temperature  at  the  instant  of  complete 
freezing  are  shown  in  Figure  II  for  e  =  0. 1  and  0.5.  For  comparison,  the  corresponding  numerical 
results  of  Tao  ( 1 967)  are  indicated  by  crosses.  Even  at  £ = 0. 1 ,  there  exists  some  discrepancy  between 
the  perturbation  and  the  numerical  solutions.  As  discussed  by  Pedroso  and  Domoto  (1973d),  and 
Stephan  and  Holzknecht  ( 1 974),  the  error  most  probably  lies  in  the  numerical  solution  itself  because 
for  e  =  0. 1 ,  one  would  expect  the  perturbation  solution  to  be  accurate. 

The  results  for  freezing  time  appear  in  Figure  1 2,  and  compare  well  with  the  numerical  predictions 
of  Tao  ( 1 967).  It  must  be  kept  in  mind  that  Tao’s  freezing  time  results  are  believed  to  be  more  accurate 
than  his  temperature  results. 
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Pritulo’s  method.  Since  its  original  exposition  in  1949,  the  spirit  of  Lighthill’s  technique  has 
remained  unchanged,  but  several  approaches  have  been  developed  to  simplify  the  procedure.  One  such 
simplification  was  introduced  by  Pritulo  (1962).  He  proved  that  the  straining  expansion  can  be 
substituted  directly  into  the  ill-behaved  perturbation  expansion  and  straining  functions  can  then  be 
chosen  in  accordance  with  Lighthill’s  condition.  In  this  way,  one  solves  algebraic  rather  than 
differential  equations  to  determine  the  straining  functions.  We  illustrate  the  procedure  by  considering 
the  inward  spherical  solidification  example  again. 

From  eq  1 83,  the  two-term  perturbation  solution  for  u  is 


1-1/r 

M  = - +e 

- 1 

1 

+ 

KJ 

1  -  1/rf 

L6(l-rfr  ' 

’’  6rf(  1  -  Tf)^ 

(237) 


To  render  eq  237  uniformly  valid,  let  us  introduce  the  expansions,  eq  1 90  and  191,  directly  into  eq  237 
and  retain  terms  up  to  0(e).  To  this  end,  we  need  expansions  for  1/r  and  1/rf.  These  can  be  obtained 
using  eq  190  and  191  and  carrying  out  the  binomial  expansions.  Thus 


i=-^--e^+o{e2) ;  J-=J — e^+ofe^) 

'■0  <1,2  V  ,,,2 

Consider  now  the  zero-order  term  in  eq  237.  It  can  be  written  as 


-  1 -!/<!»,  c  qi 
l-l/V  [<t>2(l -l/v) 


v2(l-l/\|/)2 


+  0(e2) 


For  the  first-order  term  in  eq  237,  we  simply  need  to  replace  r  by  (j,  and  r{  by  ij;.  Thus  the  two-term 
solution  for  u  is 


u  = 


1-1/<|)^  g  Oi _ ai(l  -l/<i>)  ^  \|>^-3\|/+2 

1-1/\|/  <1)^(1 -1/y)  \|/2(i  - 1/\|/)2  6(l-\j/f 


^$i^3$  +  2 
6v(l-  VF 


+  0  (e^) 


To  determine  a,  we  set  the  term  in  square  brackets  in  eq  238  to  zero,  giving 


Cl _ q|(l  -  ^^<^>)  ^  \|^-3\tf  +  2  /  j  _  i\  ^  <|>2-3(t>  +  2 

<t)2(l-l/y)  \  <])/  6\j>(l-\|/f 


=  0 


which  can  be  rearranged  as 


(238) 


(239) 
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(240) 


Oi  <l>(<l>-2)  1  \|>^-3v+2  9  1 


where  it  is  noted  that  the  right-hand  side  of  eq  240  is  a  function  of  y  alone. 
Differentiating  eq  240  with  respect  to  (|),  the  differential  equation  for  Oj  is 


o^- 


20-1  2(|.((|.-1)^ 

o(o-l)  6\|^(l-\|ff 


Integrating  eq  241  we  have 


1 

0(0-1  ) 


ai  = 


(241) 


Imposing  the  condition. 


*  *"*1  (<^1  =  0  .  makes  C  vanish.  Thus 

<|>  -^  1 


6n^(1- VF 


(242) 


which  agrees  with  eq  224. 

In  an  analogous  manner,  it  is  possible  to  derive  the  solution  for  O2  but  the  algebra  becomes  lengthy. 
The  reader  who  is  not  overwhelmed  by  the  algebra  can  verify  that  the  solution  for  03  is  in  full  agreement 
with  eq  227. 


7,2  Inward  cylindrical  solidification  with  constant  wall  temperature 

As  a  second  example,  we  consider  the  cylindrical  solidification  problem  of  section  5.2.  Following 
the  procedure  of  the  previous  section,  we  obtain  the  following  results,  which  are  taken  from  Asfar  et 
al.(1979). 


Zero-order 


A((;,^\  =  0,  Mo(v,<|>  =  l)  =  0,  Mo(v.<l>=v)=l  (243) 

d(J>  \  d<t>  / 

with  the  solution 


(244) 


First-order 

Using  the  zero-order  solution  and  making  U|  identically  zero,  the  first-order  problem  reduces  to 


9(|)  \  If  y 


(245) 
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One  condition  on  Oj  is  provided  by  eq  (192),  that  is,  We  impose  the  other  condition  by 

putting  j  ^  =  0.  The  latter  permits  the  simplest  solution  for  Oj .  The  solution  is 


^^_()>[(l-(|>^)  +  (l+(|>^)ln(|)] 
4\|/^ln^t|/ 


(246) 


Second-order 

Using  eq  244  and  246 and  making  U2  identically  zero,  the  second-order  problem,  after  considerable 
manipulation,  finally  leads  to 

— (t*— (-)!  = - ^ - ((l-«l>^)(l-9(t)^)+20(!)^ln(D 

3(1)1  3(t)\<|)//  i6(|)\|f‘*ln'‘\|f 

-36)|)'‘ln)l)  +  12)l)^ln^(i>+40(^'‘ln^<|)) - 1 - (({) +(|)ln(}) 

2y'*ln^Vf 

+  (t)^(ln(t)  -l))  + - J - (j)  In  (l>(l  -  t|;2_2\|/2(in  y  -  l)ln  y).  (247) 

2\|(‘*ln^t|r 

Imposing,  as  for  the  first-order  problem,  the  conditions  =  o  ^*^3  j  =  0,  the 

simplest  solution  for  02  is  obtained  as  I 

0^2  = - *—((t)(l-<t)^)  +  (t)(l+(t)^)ln(t))(l-v2 +  2^2(1  -lnv)lnv) 

8\|/^ln^\); 


( 3<t)  ( 1  -  (t)"*)  +  2(1)  ((t)'^  +  4<t)%  1  )l  n(l)  1 


64\|t‘*ln^V|/ 


{(21(|)^-24(t)^ 


+  3(1)) 


128y‘*lnV 

-(38(1)'^  +  8(1)^- lo)(t)  ln(l) +(20(1)'*  +  24(1)'+4)(|)  ln“(|) } . 
The  equation  for  the  freezing  time  up  to  0(e^)  is 


(248) 


A=_L 

(A|r  3^ 

3(1) 


1  +  I  (<!>.¥) 


\  <*V 


(t,=v/ 


41^  ¥ 


+  e2  /^2(¥.¥)  ,  9q2(<l>^v) 

\  <A|/  3(1) 


,  ^l(¥-¥)  ^1  (<!>■¥) 


.J1 


(249) 
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Substituting  eq  244, 246 and  248  into  eq  (249)  and  using  the  initial  condition,  y  =  1 ,  x = 0,  the  solution 
for  T  is  obtained  as 


x  =  i(l-\|f2)  +  L\i>2in\j/  +  e[li — ^ 

4  2  \  21nv  / 

+ - ^ - £^{i5(i-v2)^  +  21(  1  -^j^^)ln^J/ 

96v2ln'^\|f 

+  12(1 +\|r*)ln^\tf  +  3(l-v^)ln^\tt}  (250) 

Sample  results  for  the  temperature  distribution  atthe  instant  of  complete  freezing  (rf =0)  are  shown 
in  Figure  13.  Figure  14  shows  the  results  for  freezing  time.  For  comparison  the  corresponding  results 
of  Tao  ( 1 967)  are  also  shown.  In  Figure  1 3  even  at  e =0.01  there  exists  some  discrepancy  between  the 
perturbation  and  numerical  temperature  distributions.  However,  the  freezing  time  solutions  are  in 
good  agreement.  As  discussed  by  Pedroso  and  Domoto  (1973d)  and  Stephan  and  Holzknecht  (1974) 
the  error  in  the  temperature  distribution  lies  most  probably  in  the  numerical  solution  itself  because  for 
8 = 0.01 ,  one  would  expect  the  perturbation  solution  to  be  accurate.  Computational  experiments  with 
the  perturbation  solution  show  that  a  valid  solution  is  obtained  up  to  about  e = 0.8.  Beyond  this  value, 
the  values  of  and  (|)  obtained  gave  an  unrealistic  solution  for  u. 


Figure  13.  Temperature  distribution  for  inward  Figure  1 4.  Inward  freezing  radius  for  cylindrical 

cylindrical  solidification.  system. 


73.  Inward  spherical  solidiflcation  with  convective  cooling 

The  problem  is  described  by  eq  174  and  175  except  that  the  convective  boundary  condition  is 


where  Bi = hRJk  and  T„,  in  the  definitions  of  u,  x,  and  e,  is  to  be  interpreted  as  the  coolant  temperature 
and  not  the  wall  temperature. 

For  the  solution  of  this  problem,  Milanez  and  Boldrini  ( 1988)  followed  the  procedure  detailed  in 
section  7. 1  and  found  the  zero-order  and  the  first-order  problems  and  their  solutions  to  be  as  follows: 

Zero  order 


^(<l>“o)  -0 
a<t>2 


d<|> 


^*1 


«o(4>  =y,v)=i 


«o  = 


1 

Bi  +  (l  - 


( 1  -  Bi)\\i  +  ^ 


First  order 


1  9^(<l>ui)  _2a  2  ^uq  ^1  _2  ^'<0 

<t>  a<t>^  <t>^  <t>  9<t>  9(j)  a<i>^ 


duQ  cP^O\  _  3uo  aup 

a<i)2  ai|/ 


<l'=V 


(252) 

(253) 


(254) 


(255) 


«l(¥.  <I>=1)  =  - ’  Mi{<l)  =  y.v)  =  0 
\d<t)  d<l>  d(b  4=1 


(256) 


By  choosing  O]  such  that  =  0,  the  equations  for  Oj  and  its  boundary  conditions  are 

duQ  d^oi  ^  duQ  duQ  ^2  3^0]  2(i)^'+2a|-  ~ 

d<t»  d(t»^  d\i>  dy  4,=  y  d(|)2  d<l>  1  _  g/j2 


(257) 


a,(<t)=l.  v)=0;  ^((ti  =  l.v)=0 

do 


(258) 


and  the  solution  for  O]  is 

Oi  = - ^ - f- i (2  +  Bi)())  + J-(l  +  Bi)(l)--  -L  Bi<t>^-  i(l-B()<t)‘‘ 

t(/2[(l- tj/)B/+ B/]^  ^  ^  2  2  6- 

Note  that  as  Bi  — » «>,  eq  254  and  259  reduce  to  221  and  224,  respectively. 

The  solution  for  the  freezing  time  X  is  given  by 

x  =  Bo(^'-¥)  +  E*i(fi'.v)+o(e^) 

where 


(259) 


(260) 
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(261) 


*o(l.  v)  =  ^(l-v^);  6i(l.  v)=l-2v+ 
and  when  1, 

ho(Bi.v)=Ml-v2)+i^(l-v3)  (262) 

2  35i 


V|/)  = 


2Bi 


rV  + 


3(1  -  Bi)2  3(1  -  Bi)  9(  1  -  Bip[v(l  -  Bi)+  Bi] 


(263) 


Nowtbatthesolutionsforxandrfare  known,  the  speed  of  the  fTeezinglTont,drf/^  =  drf/d\\f(dx/d\^)~^ 
can  be  evaluated. 

Figure  1 5  shows  a  plot  of  ( 1  -  rf)  versus  x  fore = 0.5  and  for  Bi = 1 , 2, 5,<».  The  dashed  lines  represent 
the  two-term  perturbation  solution  while  the  solid  lines  give  the  numerical  solution  of  Milanez  and 
Ismail  (1984).  The  accuracy  of  the  perturbation  solution  is  excellent.  Figure  16  show  the  results  for 
e = 2  and  the  comparison  with  the  numerical  solution  indicates  that  die  perturbation  analysis  predicts 
a  slower  growth  of  the  solid  phase.  Figures  1 7  and  1 8  show  the  temperature  profiles  at  various  interface 
positions.  Once  again  the  agreement  between  the  perturbation  and  numerical  solutions  is  better  at  e 
=  0.S  than  atE=  1.0. 


Figure  15.  Freezing  time  for  inward  spherical  system  with  surface 
convection,  e  =  0.5. 


Figure  16.  Freezing  time  for  inward  spherical  system  with  surface 
convection,  e  =  2.0. 
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Figure  17.  Temperature  distribution, 
inwardsphericalsolidification  with  con¬ 
vection,  e  =  0.5,  Bi  =  oo. 


Figure  18.  Temperature  distribution,  inward  spheri¬ 
cal  solidification  with  convection,  e  =  7.0,  Bi  =  2.0. 


7.4.  Inward  spherical  and  cylindrical  freezing 
with  combined  convective  and  radiative  cooling 

Parang  et  al,  (1990)  have  considered  the  general  problem  of  inward  freezing  of  spheres  and 
cylinders  due  to  combined  convective  and  radiative  cooling.  They  have  shown  that  the  governing 
equations  for  this  problem  are 

1  +;  1  =  e  —  • — 

^  3r  \  drf  dr  drf  dr 


r=  1 ,  --i-— =  +  H(r=rf,rf)=  1  (265) 

Bit  dr 


_  I  0  (cylinder) 
‘  1  (sphere) 


(264) 


drf  _  du 
dl  dr 


r=rf 


(266) 
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where  u  =  TITf 
r=R/Ry, 
rf=R^R^ 
a  =  T^Tf 
Bi  =  hRyJk 
£  =  cT^L 
Bif  =  aeTfRyJk 
X  =  kT{t/(f>LRl) 
fi  =  BUBit  =  hl(<jeT(). 

Note  that  a  common  ambient  at  Tg  has  been  assumed  for  both  convection  and  radiation.  The  symbol 
e  represents  the  emissivity  of  the  surface.  The  perturbation  parameter  £  is  equivalent  to  the  Stefan 
number  based  on  the  absolute  freezing  temperature. 

Following  the  procedure  of  section  7. 1 ,  the  zero  and  the  first-order  problems  and  their  solutions  are 
obtained  as 

Zero-order 


±L^\+ j^=o 
3<|)  \  3<t>  /  3<t> 


(267) 


(|)  =  1 ,  _i  P(«o-a)+  (mo'*-®'*)  ;  “o(<l>  =  V,  v)=  1 
Bir  d(|) 


(268) 


For  the  sphere 
Mq  =  Cj  +  C2/<t> 

where  Cj  and  Ci  are  given  by 


C/+ 


4v 


6y^ 


(i-v)  (i-v)^ 


C,2+ 


3  I  V 

(l-Vp  Bi,(\-Vff_ 


Cl 


^  \|f>-  g^-t-  P(V-  g)-  (\|//fltr)  _Q 
(l-t|/f 

C2  =  v(l-C,) 

For  the  cylinder 
mq  =  C|  In  <|)  +  C2 
where  C\  and  C2  are  given  by 


C  4Ci^,  6Ci^  C,(4-i-p)  ^(l-g-*-i-P(l-g)  +  Ci/Bi,)_^ 
•"V  (ln\|i)^  (ln\|/)^  (iny)'* 


(269) 

(270) 

(271) 


(272) 


(273) 
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C2  =  l-C,ln\jf  (274) 

Equations  270  and  273  are  polynomials  for  the  constant  C|  and  in  general  have  to  be  solved 
numerically. 

First-order 


3<j)^  9(|> 


9^2 


9mo  duQ 
9ty  9(|) 


♦=v 


(275) 


9m  9ai  duQ 
Blr{4uQH  P) 


(276) 


Choosing  G}  such  that  U)  is  identically  zero,  the  equation  for  Gj  with  its  boundary  conditions 
becomes 


9<t)  t}f’  *y 


ai(<t>  =  l .  V)  =  0:  ^(<t)  =  l,v)  =  0 


df  ay 


The  solution  for  Gj  for  the  sphere  is 

i^(|)4+i^  <|)3_i/^+2  ^  |(|)2+ i/2^  +  3  ^](|) 
6  duf  2  chf  2\9y  dw/  6\<Atf  diit/ 


and  for  the  cylinder  is 


+  1 

2^- 

dC, 

|ln<t> 

(|,+  ^  -  ^ln(t>](t)3| 

i  <A|/  ' 

<A|/ 

<A|/ 

/  J 

V  <A|/  (A|/  <Ay  II 

(277) 


(278) 


(279) 


(280) 


T,  Time 
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Figure  19.  Freezing  radius  vs 
time  for  sphere  with  convec¬ 
tion  and  radiation. 


Sample  results  for  the  progress  of  the  freezing  front  are  shown  in  Figure  19  for  the  sphere  and  in 
Figure  20  for  the  cylinder.  The  dimensionless  ambient  temperature,  a,  is  zero  in  both  figures  and  e  is 
O.S.  The  parameter  p  is  taken  as  unity,  indicating  equal  convective  and  radiative  contributions  to 
cooling.  The  perturbation  results  match  closely  with  the  numerical  results  derived  using  the  enthalpy 
method.  However,  some  divergence  (about  1%)  is  noticeable  near  the  end  of  the  freezing  process. 
Although  not  shown  here,  the  discrepancy  between  the  perturbation  and  the  numerical  predictions 
increases  as  e  increases. 

Figure  21  shows  how  the  time  for  complete  freezing  is  affected  by  the  increase  in  the  ambient 
temperature.  This  figure  is  for  the  sphere  and  for  e = 0.5  and  P  =  1 .  As  expected,  the  total  freezing  time 
increases  as  a  increases.  Furthermore,  the  stronger  the  convective-radiative  cooling,  that  is,  the  higher 
the  value  of  Bi^  (=  Bi),  the  lesser  the  freezing  time.  As  Bi^  is  decreases,  the  difference  between  the 
perturbation  an  the  numerical  solutions  increases. 

Finally,  Figure  22  shows  the  temperature  profiles  at  the  instant  of  complete  freeze  for  e  =  0.5, 1 .0, 
and  2.0.  These  curves  pertain  to  spherical  freezing  and  correspond  to  a  =  0,  and  Bi  =  Bi^  =  1.  The 
maximum  error  between  the  perturbation  and  numerical  solutions,  which  occurs  in  the  vicinity  of  the 
center,  is  about  8%. 


Figure  20.  Freezing  radius 
vs  time  for  cylinder  with 
convection  and  radiation. 


Figure  21.  Final  freezing 
time  for  sphere  with  con¬ 
vection  and  radiation. 
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Figure  22.  Temperature  profiles 
for  sphere  with  convection  and 
radiation. 


8.  METHOD  OF  MATCHED  ASYMPTOTIC  EXPANSIONS 

The  most  widely  used  singular  perturbation  technique  is  the  method  of  matched  asymptotic 
expansions.  In  this  method  the  outer  expansion,  which  is  the  regular  perturbation  expansion,  is 
supplemented  by  an  inner  expansion  in  terms  of  a  stretched  variable,  which  then  describes  the  behavior 
of  the  function  in  the  region  where  the  outer  expansion  breaks  down.  Finally  the  outer  and  the  inner 
expansions  are  matched  to  derive  a  uniformly  valid  solution.  A  criterion  which  is  commonly  used  is 
Van  Dyke’s  matching  principle.  The  method  will  be  illustrated  with  the  help  of  two  examples. 

8.1  Vaporization  of  a  semi-infinite  solid  due  to  constant  heat  flux 

The  use  of  high  power  lasers  and  electron  beams  to  induce  evaporation  (sublimation)  of  solid 
materials  has  become  important  in  many  materials  processing  applications.  By  ignoring  the  interme¬ 
diate  melting,  the  complexity  of  the  problem  is  considerably  reduced. 

Consider  a  planar  solid  of  semi-infinite  extent  heated  by  a  constant  heat  flux  q,  applied  at  x = 0.  In 
the  simplest  case  one  can  envision  that  all  the  energy  absorbed  at  the  surface  is  used  to  vaporize  the 
material  and  none  is  conducted  into  the  solid.  This  vaporization-controlled  limit  can  arise  if  the  heat 
is  applied  too  rapidly  for  conduction  to  occur.  If  the  temperature  distribution  ahead  of  the  vaporizing 
boundary  approaches  a  steady  state,  the  velocity  will  also  attain  a  constant  value.  In  this  situation,  the 
velocity  of  the  vaporizing  front  v  can  be  found  as 


pWTv-t’O-hl] 

where  the  extraction  of  the  sensible  heat  has  been  accounted  for. 

However,  a  more  realistic  model  must  include  the  effect  of  transient  heat  conduction  into  the  solid 
ahead  of  the  vaporizing  boundary.  It  is  this  model  that  is  developed  here  and  analyzed  using  the  method 
of  matched  asymptotic  expansions.  The  one-dimensional  transient  heat  conduction  equation  and  the 
associated  conditions  are 

d^T  _  1  dT 
dx^  ®  dt 
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(282) 


r(xy,r)=r^;  =pL^-q  (283) 

=5  JCv 

T{x,0)=Ti-,  T{oo,t)=Ti 

where  Xy  denotes  the  location  of  the  vaporizing  boundary  and  Ty  is  the  vaporization  temperature  of  the 
solid.  It  is  assumed  that  the  solid  is  initially  at  temperature  T^,  and  that  at  large  distances,  the 
temperature  of  the  solid  remains  unaffected. 

Introducing  the  following  dimensionless  quantities: 

0  =  (r-7’i)/(rv-7’i),  X  =  xv/a,  \  =  xyv/(x  (284) 

x=v^r/a,  e  =  c(Ty—Ti)/L 

into  eq  282  and  283  gives 


_  89 
dX^  8t 


(285) 


0{A„t)=l  ;  f^-1  l-e(^  +lUo  (286) 

\dx  }  \dX  x=\  / 

0(X,O)  =  O;  0(oo,x)  =  O 


Let  us  assume  a  regular  perturbation  solution  for  0,  X,  and  T)  =  ^  as  follows: 

dl 


0=  00+  e  01  +0(e2) 

(287) 

X=  Xo+  eXi  +0(e2) 

(288) 

Ti  =  qo+  eT|i+0(£^) 

(289) 

Substituting  eq  287-289  into  eq  285  and  286,  we  obtain  the  following  equations  for  0^  and  T)o: 

II 

O 

(290) 

^0=1  or  ^=1 
</t 

(291) 

X=Xo,  00=1  . 

(292) 

Ignoring  the  preheating  effects  while  the  boundary  is  being  raised  to  its  evaporation  temperature, 
since  such  effects  are  important  only  for  x = (K£^)  as  seen  later,  we  can  integrate  eq  29 1  with  the  initial 
condition  T  =  0,  Xq  =  0  to  give 

II 

(293) 
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To  solve  eq  290,  we  first  make  the  transformation  from  die  fixed  cowdinate  system  (X,t)  to  the 
moving  coordinate  system  (X,  x)  and  rewrite  eq  290  in  moving  coordinates  (see  Carslaw  and  Jaeger 
1959,  p.  13) 


d^§o  +  3Qo  _  9Qo 
9X  9x 


(294) 


Taking  the  Lqilace  transform  of  eq  294  and  noting  that  6(X,0)  =  0,  we  get 

+  0o(j)  =  o  (295) 

dk 

where  6o(s)  is  the  Laplace  transform  of  6  and  s  is  the  Laplace  variable.  The  boundary  conditions  for 
eq  295  are 


X  =  0,6o=l  or  X  =  oo,  e^j=:0or  0o(it)  =  O  . 

The  solution  of  eq  295  subject  to  eq  296  is 


(296) 


S 

From  Carslaw  and  Jaeger  (1959),  p.  495,  the  inverse  of  eq  297  is  obtained  as 

^  (X  ,x) = i  erfc  ^  ^ 

2  \2Vx/  2  \2Vx/ 

Since  the  solution  for  X  to  (Ke)  is  given  by  eq  293,  we  can  write  X  =  X  -  x,  giving 


(297) 


(298) 


We  now  consider  the  first-order  problem  for  X,  which  follows  as 


(299) 


+  1 
x=k) 

Differentiating  eq  299,  we  have 


^^_dXi  _89o 
3x  9X 


(300) 


IX-x 

-  erfc  2 - .  J-  e 

Vx  2 


(301) 
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Further  noting  that  erfc  (-X)  =  2  -  erfc  X,  eq  301a  can  be  expressed  as 


dX 


=-l  + 


x=V) 


J-erfc|-Lx*^|-(nx) 


(302) 


Using  eq  302  in  300,  the  two-term  solution  for  q  becomes 


q  =  l  +e 


-i-  erfc  X  j  -  ( n  x)“  c 


(303) 


Equation303  isagoodapproximationforx>  1  butforx=0(e^),theterm-e(nx)~’^  e  “^becomes 
(K 1 )  rather  than  (Ke)  which  is  an  indication  of  eq  303  breaking  down  when  x = 0(c^),  diat  is,  for  small 
times.  Thus  the  solution,  eq  303,  constitutes  an  outer  solution  valid  for  long  times.  For  short  times,  an 
inner  solution  must  be  constructed. 

Inner  Solution 

The  preheating  effect,  which  was  ignored  in  developing  the  outer  solution,  becomes  important  in 
determining  the  motion  of  the  vrqrorizing  boundary  for  small  times,  x  =  0(e^).  During  the  preheating 
time,  the  heat  conduction  (eq  285)  is  subject  to  die  following  boundary  and  initial  conditions: 


dx 


=-l-  -I-;  0{<»,  x)  =  0 


(304) 


lx=0 

0(X,O)  =  O 

where  the  first  condition  in  eq  304  is  obtained  from  the  constant  heat  flux  condition  at  X = 0,  that  is, 

3r| 


-k 


dX 


=  <!■ 


x=o 


The  solution  to  this  problem  can  be  adapted  from  Carslaw  and  Jaeger  (1959)  and  expressed  as 

(305) 

The  preheating  time  can  be  obtained  by  using  0  - 1  and  X  =  0  and  is 


e=(i  +i)[2(T/,r  «p(-^)-Ar«fc(^)]. 


Xp  =  i  Jt£2/(  1  +  e)^  (306) 

Equation  306  shows  thatXp =0(e^)  and  substantiates  our  earlier  assumption  about  the  preheating  time. 
If  we  substitute  eq  306  into  eq  305,  we  obtain  the  temperature  distribution  in  the  solid  at  the  end  of  the 
preheating  period  as 


e(X,Tp)  =  exp[-{x(l+  e)/(n'^e)}2]-{X(l+  e)/e}erfc{x(l+  e)/(jt'/2e)}  (307) 

Examination  of  eq  306 and  307  shows  that  the  appropriate  stretching  of  the  variables  for  small  times 
is 


X*  =  X/E,  =  Tp)/e2,  =  e*=0,  11*  =  T1 


(308) 


Substituting  the  transformation,  eq  308,  into  eq  285  and  286,  the  equations  for  the  inner  region  become 

d^e*  ao* 


ax 


*2 


dz* 


(309) 


AtX*  =  ^=eX*,  0*  =  1,  ■n*  =  l+-^^+e 

e  ax* 

At  x*  =  oo,  e*=o 

At  x*=0,  0*  =exp[-{x*(l  +e)/n*^)  ]-{x*(l  +e) }  erfc  | . 

1  I 

Again  we  seek  a  perturbation  solution  of  eq  309-312  in  the  form 

0*=0S+  eOj 

tl*  =  Ti5+  eti* 

The  zero-order  problem  is  governed  by 

,  30^ 
ax  *2  ax* 


(310) 

(311) 

(312) 

(313) 

(314) 

(315) 


AtX*=0,  05=1  iiS=l+M  (316) 

ax 

x*  =  oo,  e5=o  (317) 

x*  =  0,  0o  =  exp(-X*^/7t)- X*erfc(X*/i^)  (318) 

Since  the  times  involved  are  0(e^),  the  motion  of  the  vaporizing  boundary  during  this  time  would 
be  very  slight.  Thus  one  can  ignore  the  last  condition  in  eq  316  and  treat  the  problem  as  a  no-phase- 
change  heat  conduction  problem.  Taking  the  Laplace  transform  of  eq  315  we  get 

^-^^-j05(j)  +  exp(-X*^/7i)-X*  erfc(x*/V7t)  =  0  (319) 

ax*2 
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where  0o(5)  is  the  Laplace  transform  of  0o  and  s  is  the  Laplace  variable.  The  boundary  conditions,  eq 
3 1 6  and  317,  transform  to 


X*  =  0,  eS(j)=i;  X*==o,  0;(s)  =  O 


(320) 


The  solution  of  eq  3 1 9  and  320  using  the  variation  of  parameters  method  can  be  obtained  as 

ii/2 


0o(j)  =  J-exp  2exp(-s'^X*)  erf|j-5nj* 

+  exp(5*^X*)erfc  |x*/ji'^+ 


-exp(-j'^X*)erfc|x*/Jt'^- 

+  iexp(-X  ^/Jt)-(x*/j) erfc(x*/it*^) 

The  Laplace  transform  of  the  last  condition  in  eq  316  is 


-m 


(321) 


dX 


(322) 


X  =0 


Using  eq  321  to  evaluate  and  substituting  in  eq  322,  the  solution  for  t1o(j)  is  given  by 
1)5(5)=  J-  exp  |in5|erfcj^i(5ic)''^ 

The  inverse  of  eq  323  is  given  by 


(323) 


q5=—  arcsin 
K 


(324) 


Matching 

Expressing  the  outer  solution,  eq  303,  in  terms  of  the  inner  variable  T*,  and  expanding  it  to  two  terms 
and  keeping  t*  constant,  we  find  that  this  expansion  matches  with  the  one  obtained  by  expressing  the 
inner  solution  eq  324  in  terms  of  the  outer  variable  t  and  expanding  it  in  two  terms  keeping  t  constant. 
Thus  the  matching  condition  is  automatically  satisfied. 

Following  van  Dyke  (1975)  we  may  now  construct  a  uniformly  valid  solution  by  adding  the  outer 
expansion  and  the  inner  expansion  and  subtracting  the  solution  in  the  domain  of  overlap.  Thus 

q=l+  e  |J-erfc|J-x’^j- 


+  ^  arcsin 
Jt 


1  +  -^ 


4(t-  rpl 


-1/2 


1- 


(jcr) 


1/2 


TS  To 


(325) 


Equation  325  is  a  uniformly  valid  solution  having  an  error  of  0(e)  when  x  =  O(e^)  and  of  0(e^)when 
X  =  0(  1 ).  Table  6  and  Figure  23  compare  the  perturbation  solutions  to  a  heat  balance  integral  solution 
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Table  6.  Comparison  perturbation,  heat  balance  integral,  and 
numerical  solution,  constant  wall  heat  flux,  e  =  0J2257. 


Heal  balance 


a 

Perturbation 

Numerical 

integral 

Inner 

Outer 

Uniformly 

valid 

Landau 

(1950) 

Goodman 

(1958) 

0.0269 

O.OS183 

0.3316 

0.1090 

0.0805 

0.01423 

0.02782 

0.1076 

0.3441 

0.2162 

0.1768 

0.0584 

0.03634 

0.2913 

0.4388 

0.3981 

0.4015 

0.3178 

0.06S6 

0.4959 

0.6075 

0.6005 

0.6288 

0.6129 

0.20 

0.7149 

0.8135 

0.8132 

0.8357 

0.8582 

1.0 

0.8727 

0.9549 

0.9549 

0.9564 

0.9845 

9.0 

0.9576 

0.9994 

0.9995 

0.9995 

— 

100.0 

0.9873 

1.000 

1.000 

1.000 

— 

and  a  numerical  solution.  The  numerical  results  of  Landau  ( 1 950)  for  e  =  0.225  are  shown  as  circles 
for  comparison.  Figure  24  shows  a  plot  of  t)  versus  t  obtained  from  eq  325  for  e  =  0  (vaporization- 
controlled  limit),  0.05  (graphite),  0. 1 5  (tungsten)  and 0.25  (lead).  The  horizontal  intercepts  correspond 
to  die  preheating  time  Tp_  eq  306.  The  integration  (trapezoidal  rule)  of  these  curves  gives  the  location 
of  the  vaporizing  boundary  as  a  function  of  tand  is  shown  in  Figure  25.  The  curve  marked  e = 0  is  the 
limiar  relationship  represented  by  eq  293.  The  agreement  with  the  present  results  is  good,  and  confirms 
the  validity  of  eq  325. 


SJ,  Vaporization  (sublimatkm)  of  a  finite  solid  due  to  constant  heat  flux 

Here  we  consider  the  problem  of  section  8. 1  but  now  the  solid  has  a  finite  thickness  t  as  shown  in 
Figure  26.  The  face  at  x  =  0  is  assumed  to  be  insulated.  Let  us  first  consider  the  preheating  when  the 
surface  temperature  is  elevated  from  Tj  to  Ty.  For  this  no-phase-change  problem  we  have 


_\dT 
dx^  ~  «  at 


(326) 
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Figure  24.  Melt  velocity  vs.  time, 
matched  asymptotic  solution. 


Figure  25.  Melt  position  vs.  time, 
matched  asymptotic  solution. 


q 


Vapor 


x~t 

X  »  (vaporizing 
boundary) 


Solid 


Insulated 


J 


x-0 


Figure  26.  Vaporization  of  finite  solid  due  to  constant  surface 
heat  flux. 
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x  =  0,  — =0;  x  =  t,  k—=q;  t=0,T=T,. 

dx  dx 


(327) 


The  solution  of  eq  326  and  327  is  given  by  Carslaw  and  Jaeger  (1959)  as 


0  =  j2x  +  2  — — i^lii-cos(>iiiAr)exp(-n2n^x) 


(328) 


where  0  =  (T-Tij/fTy-Ti),  jc  =  Jt/f.  x  =  at/f^.and  g=— — ^ The  preheating  time  Xp  when 

A:(Tv-Ti) 

6  =  1 ,  X  =  1 ,  is  given  by  the  following  implicit  relationship 


Xp  + 


i--i  £  -L  exp(-n27t2xp) 


3  n2„t'in2 


-L 

Q 


(329) 


Next,  consider  the  phase  change  problem.  The  heat  conduction  into  the  solid  is  still  governed  by 
'•q  326  but  the  conditions  to  be  satisfied  are 


jr=0. 


—  =  0;  jr  =  jrv.  T=  Ty  ; 
dx 


X 


Xyj  , 


(330) 


The  initial  condition  is  given  by  eq  328  with  x  =  Xp. 

Introducing  the  following  dimensionless  variables  :(► = (Ty  -  T)l{Ty,  ~  Tj),  q  =  x/xy,  <y=Xy/i,  e = ciTy 
-  Ti)/L,  2  =  £x  =  eou/^2  into  eq  326  and  330,  we  obtain 


^^.  =  e[o2^-qo^.  ^ 
dz  dq. 

(331) 

q  =  l,  <|>=0;  q=0,  ^  =  0 
dq 

(332) 

z  =  0,  (|>=y(q) 

where/Iq)  is  given  by 


Aq)=l-GXp-G 


cos(n7tq)  exp(-q2jt2xp) 


noting  that  x  =  Xp,  Xy  =  f,  a  =  1,  q  =  X,  and  <|>  =  1  -  0. 
The  "nergy  balance  at  the  interface  transfonns  to 


(2<y+— ( I  ,x)=-J-a^=-o^ 
dq  E  dx  dz 


(333) 


(334) 
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Outer  solution 

The  long  time  or  the  outer  solution  is  assumed  as 

<|>  =  (1)0  +  e)t)i  (335) 

0  =  00  +  £<^1  (336) 

Substituting  eq  335  into  eq  331  and  332  and  equating  coefficients  of  e^,  we  have 

^^  =  0  (337) 

dn" 

q=0,  ^  =  0;ti  =  1,  (t)o  =  0  (338) 

dn 

The  solution  of  eq  337  and  338  is 

(|)o  =  0  (339) 

Equation  339  shows  that  after  a  long  time,  the  solid  is  at  a  uniform  temperature  throughout.  Since 
(jio  =  0,  it  is  easy  to  see  that  ([i,  =  <1)2  =  ...(l)n  =  0. 

Substituting  eq  336  into  eq  334  and  equating  coefficients  of  e®  the  equation  for  Oq  is  obtained  as 

^  =  (340) 

dz 

Integrating  eq  340  we  have 

Oo  =  -Cz  +  C,  (341) 

where  the  constant  C|  would  be  determined  by  matching  the  outer  and  the  inner  solutions. 
Equating  coefficients  of  e,  gives  the  following  equation  for  O] : 

ao^  +  o,^  =  -ea,  (342) 

dz  dz 

Using  eq  341  for  Oq  in  eq  342,  it  can  be  seen  that 

^  =  0  or  ai=C2  (343) 

dz. 

where  C2  is  a  constant  which  would  be  determined  by  matching. 

Inner  solution 

Examination  of  eq  331  reveals  that  for  the  inner  solution,  the  appropriate  time  variable  is  x  =  s/e. 
Thus  the  inner  expansion  is  written  as 

<t)*  =  <|)S+ e<t)i  (344) 

o*  =  o5+ea*  (345) 
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where  (t>*  =  (|>i  and  o*  s  a. 

Substituting  eq  344  and  345  into  eq  331  and  332  and  equating  coefficients  of  e^,  we  obtain 

•2  ^  3^ 

3x  (A  3q 

q  =  l.  <t>o=0,  ^=0;  T  =  0.  <\>o=m) 

dn 

Similarly,  substituting  eq  345  into  334  and  equating  coefficients  of  e~',  we  get 

ao^=0  or  ^  =  0  (348) 

dr  dr 

Since  at  x  =  0,  a  =  1 ,  the  initial  condition  on  Oq  is  x  =  0,  Oo  =  1 .  Thus  the  solution  of  eq  348  is 

aS=l  (349) 

In  view  of  eq  348  and  349,  eq  346  reduces  to 


(346) 

(347) 


3^<|>o  _  3<l>o 
3x 


(350) 


The  solution  of  eq  350  subject  to  eq  347  can  be  obtained  using  the  method  of  separation  of  variables 


as 


<t>o=  X  ^mCOs(X,mTl)exp(-^Si't) 


(351) 


m  =  0 


where  A,, 


,=  (»*!)» 


and 


{-O'" 


1  +  0j-L-Xp-i 


+  2G  y  (jlL):  Uxpl-n^Tt^Xp) 

Jt^  n=I  I  rtlt-A-m  nK+Xm  I 


To  obtain  the  equation  for  a*,  we  substitute  eq  344  and  345  into  eq  334  and  equate  coefficients  of 


e”  to  give 


dti 


■do, 


^0' 


\  dt  dr 
which,  in  view  of  eq  348  and  349,  reduces  to 


(352) 
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(353) 


dn 


Using  eq  35 1  to  evaluate  ( 1 ,  T)and  substituting  the  result  into  eq  353,  the  solution  for  at,  subject 
3t1 

to  the  initial  condition,  t  =  0,  ai  =  0,  is  finally  obtained  as 


cs\  =  -Qx  +  X  5-^/lra[l-exp(-A,^T)] 


B=0 


(354) 


Matching 

Writing  the  outer  solution  a = -  Qz + Uj  +  e  C2  in  terms  of  the  inner  variable  x  and  matching  it  with 
the  inner  solution 


o*  =  1  +  e 


expressed  in  terms  of  the  outer  variable  z,  we  find  that 


C,  =  l  and  C2=  I  (-irA^/X^. 

m=0 

The  outer  and  the  inner  expansions  agree  in  the  intermediate  time  zones.  Eliminating  the  common 
parts  of  the  two  expansions,  a  uniformly  valid  solution  is  obtained  as 


<l>=  S  '4inCOs(^m^)®^p(-^m‘f)  +  0(e)  (355) 

/n  =  0 

a  =  l+ejx  — -4n,[l-exp(-X^T)]-(?t|+0(e2)  (356) 

Equation  355  shows  that  for  long  times,  <|)  approaches  zero,  which  in  turn  implies  that  T approaches 
Ty.  For  large  times,  eq  356  shows  that  o  -» ( 1  -Qx).  For  complete  vaporization,  a  =  0  and  assuming 
that  at  that  time,  the  exponential  term  becomes  negligible,  the  total  time  for  complete  vaporization,  Ty 
is  given  by 


1+e  X 

ni=0 


(357) 


where  the  second  term  in  eq  357  represents  the  effect  of  transient  heat  conduction  on  the  vaporization 
time.  The  first  term,  l/(eQ),  represents  the  vaporization-controlled  limit.  The  effect  of  transient 
conduction  is  to  increase  the  time  tocomplete  vaporization  because  some  of  the  energy  supplied  is  used 
to  heat  the  solid  rather  than  being  used  for  the  vaporization  of  the  solid. 
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9.  EXTENDED  PERTURBATION  SERIES  METHOD 


The  perturbation  expansions  that  we  have  discussed  so  far  were  al  I  termi  nated  at  the  second  or  third 
term.  Such  a  truncated  expansion  is  valid  for  a  limited  range  of  values  of  e.  If  used  beyond  the  range 
of  applicability,  the  approximation  fails  to  converge  and  gives  erroneous  answers.  The  method  of 
extended  perturbation  series  attempts  to  remove  this  limitation  by  combining  analysis  and  numerical 
computation. 

The  method  follows  three  steps.  First,  the  set  of  perturbation  equations  is  programmed  for  solution 
on  a  digital  computer  so  that  a  large  number  of  temis  can  be  generated.  Second,  the  coefficients  of  the 
series  are  utilized  to  identify  the  location  and  nature  of  singularities  limiting  the  range  of  applicability 
of  the  series.  With  this  knowledge,  the  final  step  is  to  recast  the  series  using  one  or  a  combination  of 
devices  such  as  the  Euler  transformation.  Shanks  transformation,  Padd  approximants,  extraction  of 
singularities,  and  series  reversion.  The  improved  series  generally  has  better  accuracy  and  a  wider  range 
of  applicability  than  the  original  series. 

9.1  Extensioii  of  series 

For  simple  problems,  the  addition  of  a  few  extra  terms  to  a  two-  or  three-term  expansion  can  be 
achieved  with  hand  computation.  However,  complex  problems  require  the  use  of  a  digital  computer. 
Depending  upon  the  problem,  there  are  two  possible  approaches.  First,  if  the  problem  is  such  that  a 
pattern  can  be  established  for  the  sequence  of  solutions,  then  one  can  write  a  program  incorporating 
this  pattern  and  calculating  the  terms  in  sequence.  However,  to  establish  the  solution  pattern  it  is 
essential  to  calculate  the  first  few  terms  by  hand.  These  hand  computations  also  assist  in  debugging 
the  program.  Second,  if  the  solution  pattern  is  not  discernible,  one  must  adopt  a  fully  numerical 
procedure  to  solve  the  sequence  of  perturbation  equations. 

An  example  of  the  former  approach  in  phase  change  heat  transfer  is  provided  by  Pedroso  and 
Domoto  ( 1 973a),  who  give  the  details  of  how  the  terms  of  the  perturbation  series  for  eq  42  and  43  can 
be  generated  automatically.  Using  their  nomenclature,  these  equations  are  vmtten  as 


dx'^  9-rf 


(358) 


M(x=Xf,A:f)=l,  H(A:=0,A:f)=^  ,  ^  (359) 

x=Q  ^  x=xi 

where  u  =  {T-  Tq),  x  =  hX/k,  Xf = hXffk,  x = h\Tf  -  TQ)xJpLk,  e  =  c{Jf-  Tq)IL  and  X  denotes 

the  dimensional  distance. 

The  solution  of  eq  358  is  assumed  to  have  the  form 


m(x,  Xf;  £)=  Mi(x,  A:f)e'~',  i  =  1,  2,  ••  •, /V,  (360) 

where  the  summation  convention  over  repeated  indices  is  used  in  eq  360  and  is  the  number  of  terms 
in  the  perturbation  series.  Substituting  eq  360  into  eq  358  and  359  and  equating  coelficients  of  equal 
powers  of  e,  one  obtains 


^=0  for,  =  l 
dx^ 
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_  ^ 

"  dx 


X=Xf 


j  =  \x  ••■.i-l  for  iS2 

OtCf 


Ui{x  =  Xi,Xf)  = 


1  for  i  =  1 
0  for  i  S  2 


ox 


i=0 


(361) 


The  solution  to  the  differential  equations  in  eq  361  can  be  obtained  as 


(362) 

where  (  =  1,2,...  Ni,j=  1, 2, ...  2(,  andXjj  j  are  functions  of  JCf  to  be  calculated  from  the  boundary 
conditions  in  eq  36 1 .  Pedroso  and  Domoto  ( 1 973a)  show  that  ]  can  be  evaluated  from  the  following 
relationships. 

^1,1,1  =  —  (363a) 

1  +  Xf 

A.U4n  =  >-134n  =  (-l)”~'  (363b) 

(l  +^f)'" 

where 

m  =  2,3,. no  sum  on  m. 

(363c) 

/(/+l) 


where 

i  =  2,3,...,lV,;Z=l,2,...,2(i-l) 


/=1A 


2 


*=2,3, 


•,2( 


no  sum  on  1. 


X,u  =  A.ai=^^ 

1  +  X{ 


i-l 


where 


(363d) 


(=2,3,...,ZV,;  7=3,4,. ..,2/ 


^i,/+2,m  ~  ^ni,m-n+l 


— - - 


(363e) 


where 


1=2,3,  /  =  1,2,.  .,2((-I) 
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2 


ii+n 


.  k=23.  -  -.2i 


m=2^,  i  +  \‘,  «  =  1A  •••.»» 


A  =  I^.-.,n,  nosumonin 

^Ijn”  ^,7,m  —  ~  Anjn-n  +  l  <4m-h  +1  ^ijj».^-b+l  ^Ujn-a  +I 
whoe 

I  =  23,  -,A4,-1  ;  y=3.4.  --,2/ 
m  =23.  •  *  +  1  ;  n=13,---,»n 

A  =  13.---.n.  nosiunon/iL 

The  quantities  n  and/are  defined  as  follows 

=  ^^j  =  1  fw  fsl.  3.-.-  nosumoni 
qj  =  <3i_I.j_I  +  q_t,  jfor/=3,4,  >  =  23, 

4i=-*r‘ 

and  the  (m-l)th  derivative  of  the  function /is  given  by 

(0- 1 )  0-2)-  •  ■  (/« -  m  +  Djc/"'"  for  m  =  2.3,. .  \ 

/ 

OformSj  +  1  nosumonyorm.  f 

Combining  eq  359c  and  360.  the  equation  fortfc/drf  can  be  obtained  as 

dxf  dxf 


/-Lfori  =1 


-  ^  for  f  =  2 

--L(gi+  gj^,  ^)for.  =  3.4.  ...N, 
dtf  / 


|whcrey=  13.  ■■■i-2 


& 


where  i  =  1,2,...  Al, ,  >  =  2,3 . 2i . 


(363f) 


(364) 

(365) 


(366) 


(367) 
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For  a  given  value  ofxf,  eq  367  can  be  evaluated  once  the  values  of  |  are  known.  By  integrating 
the  terms  dx  jdxf ,  one  can  obtain  the  coefficients  tj.  In  section  3.2,  the  coefficients  Xj  were  denoted 
by  c„  and  Table  1  contains  the  first  nine  values  for  a  range  of  Xf  values. 

An  example  of  the  extension  of  the  perturbation  series  using  a  fully  numerical  procedure  is  given 
by  Beckett  (1981).  He  considered  the  inward  freezing  of  a  cylinder  and  extended  the  series  for  the 
freezing  front  location  E,  in  time  e,  to  31  terms.  For  p  =  0.1  (Stefan  number  =  10)  his  result  is 

£■=2.5139  +0.5130  e +  0.6540  £3^  +  1.2779  £2  +  2.2475 

+  4.8757  £3+11.183  +  26.611  £'^  + 65.528  £’^+164.64  £5 

(368) 

+  421.22  £"'2  +  1093.6  £6 +  2874.5  £’3^2  +  7643.2  £2+... 


+  8.651  X  10^£21'2  + . . .  +  2.759  x  10“  £3'/2 . 


9J1  Analysis  of  series 

Once  the  extension  of  the  series  has 
been  accomplished,  the  next  step  is  to 
explore  the  analytic  structure  of  the 
solution  identifying  the  location  and 
nature  of  singularities,  if  any  exist. 
The  final  pattern  of  the  sign  prevailing 
in  the  series  determines  the  location  of 
the  dominant  singularity.  If  the  signs 
are  fixed  as  in  eq  368,  the  singularity 
lies  on  the  positive  axis,  but  if  the  signs 
alternate,  it  lies  on  the  negative  axis. 
With  random  signs,  the  most  likely 
possibility  is  that  singularities  occur  as 
complex  conjugate  pairs. 

When  the  nearest  singularity  lies 
on  the  negative  axis,  it  usually  carries 
no  physical  significance.  On  the  other 
hand,  a  singularity  appearing  on  the 
positive  axis  can  often  be  interpreted 


plot. 


physically.  For  example,  the  positive  axis  singularity  associated  with  the  series  in  eq  368  indicates  the 
completion  of  the  freezing  process.  In  some  cases,  the  positive  axis  singularity  points  to  the  limit  of 
validity  of  the  mathematical  model  itself.  Another  possibility  with  the  positive  axis  singularity  is  that 
it  is  not  real  but  simply  an  indication  that  the  function  is  double  valued. 

To  establish  the  location  and  nature  of  singularities,  the  best  approach  is  to  calculate  the  ratio 

Icn/Cn_il  for  the  series  ^  c„£'' and  plot  it  against  l/n,  giving  what  is  called  a  Domb-Sykes  plot.  For 
n  =  0 

large  n,  the  ratio  Cn/c„.]  often  becomes  linear  in  I/n  and  is  of  the  form 


^n-l 


(369) 


where  ( 1  +  a)  is  the  slope  and  £o  is  the  radius  of  convergence  of  the  series.  By  extrapolating  the  graph 
to  l/n  =  0,  one  can  obtain  the  intercept  I/Eq.  Figure  27  shows  the  Domb-Sykes  plot  for  the  series  eq 
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368  from  which  a  =  3/2  and  Eg"*  =  8.77.  Thus  one  can  estimate  the  dimensionless  time  for  complete 
freezing  as  ep  =  0. 1 1 40. 

93  Improvement  of  series 

As  the  final  step,  our  knowledge  about  the  leading  singularity  is  used  to  improve  the  series.  A 
number  of  improvement  devices  are  available.  The  choice  of  a  particular  method  depends  on  the 
direction,  distance,  and  nature  of  the  singularity  as  revealed  by  the  Domb-Sykes  plot.  Often  the 
problem  is  such  that  more  than  one  device  is  applicable,  and  it  is  always  enlightening  to  try  different 
possibilities.  Also,  if  a  single  technique  is  not  adequate,  one  should  consider  using  two  or  more  of  them 
in  combination. 

Euler  transformation.  When  a  perturbation  series  contains  alternating  signs,  the  singularity  lies  on 
the  negative  axis  and  carries  no  physical  significance.  In  this  case,  the  simplest  device  to  use  is  an  Euler 
transformation  based  on  the  estimate  of  Eg.  Widi  this  transformation  the  singularity  is  mapped  away 
to  infinity.  The  advantage  of  this  device  is  that  the  exact  nature  of  the  singularity  need  not  be  known. 

Let  the  perturbation  series  be 


I  £"^4. 

n=0 


(370) 


and  the  nearest  singularity  be  located  at  E  =  Eg  (estimated  from  a  I>omb-Sykes  plot).  The  transforma¬ 
tion  envisages  using  a  new  variable  e*  such  that 


e*  =  _S— 

E+  Eg 

which  gives  the  Eulerized  series  as 


(371) 


nssO 


(372) 


where  the  coefficients  b„  are 


bg  =  ao 


^  =  I 


(«-l)! 


5(n-7)!0-l)! 


(373) 

(374) 


Although  the  Euler  transformation  eq  37 1  has  been  written  for  the  power  series  in  e,  one  can  also 
use  it  for  a  power  series  in  a  variable.  Consider,  for  example,  the  series  in  eq  183  which,  for  inward 
freezing,  breaks  down  as  rf  -» 0.  Rewrite  it  as 


u  _ 
MO 


where 


ai  =  J- 

6rf 


(375) 


(376a) 
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03 


36rf  \  120  L  \rfl 


(376b) 


The  series  in  eq  375  can  now  be  regarded  as  a  series  in  1/rf.  Note  that  as  rf  0,  aiUg  and  a-^UQ  remain 

finite.  For  convenience,  we  can  express  eq  375  as 


Ji-=  1  +  Ciq+  C3q3 
«0 

.^c2. 


(377) 


where  Cj  =  eai  and  C3  =  are  functions  of  r,  r^,  and  e.  Applying  the  Euler  transformation 

q*  = 

q+AT 


(378) 


to  eq  377,  we  have  the  following  Eulerized  series: 

-U-=l  +  Kai  e(l  +  q*)q*+  ATeCoi  + 

«0 

The  series  in  eq  379  remains  valid  as  rf  — » 0  or  q  — »  00  because  as  q  ->  00,  q*  — »  1 . 

If  one  constructs  a  regular  perturbation  series  for  then  it  is  found  that 

dx 


(379) 


^=1+ eZ>i  i+ e263-I-  (380) 

80  ff  q3 

where 

f,3=i±i(t  (381) 

3  45 

If  the  transformation  eq  378  is  now  applied  to  eq  380,  we  obtain 

1  +  Kbi  efl  +  q*)q*  +  Ke{bi  +  K^bj  e)q*3  (382) 

80 

If  eq  379  and  382  are  used  in  the  overall  energy  balance  for  the  entire  duration  of  the  freezing 
process,  one  can  obtain  an  equation  for  K  =  K(£).  It  has  been  shown  by  Pedroso  and  Domoto  (1973c) 
that  if  only  terms  in  the  first  power  of  q*  are  retained  in  eq  379  and  382,  then  the  overall  energy  balance 
gives  the  following  transcendental  equation  for  K: 

+  =  (383) 

6K  \  3-eKl  40  (3-ef:) 

The  solution  of  eq  383  for  a  range  of  values  of  e  is  given  below: 

£  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  2  3 

K  14.96  9.099  6.699  5.350  4.474  3.854  3.391  3.030  2.791  2.503  1.351  0.9287 

Integrating  eq  382  to  obtain  the  freezing  time  as  a  function  of  r^,  the  solution  becomes 
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(384) 


T*  (Euler )  = -£- [3  +  ( 3  -  e)jf] 
9K 


,  3-eK ,  3  +  (3-e)A: 

1-rf - In — - - ^ 

3K  3+(3rf~t)K 


3) 

6  3 

The  results  from  eq  384  are  compared  with  the  numerical  results  of  Tao  ( 1 967)  in  Figure  28  for  e  = 
0.1, 0.5  and  1.0.  While  the  regular  series  diverged  as  rf->0  (see  Fig.  10),  theEulerizedsoiesineq 


384  agrees  quite  well  die  numerical  results. 

Shank  tran:^omuUions.  Shanks  (I9SS) 
introduced  a  family  of  four  nonlinear  trans¬ 
formations  to  accelerate  the  convergence  of 
slowlyconvergentanddivergentseries.The  t 
merit  of  these  transformations  is  that  they 
do  not  require  any  information  about  the 
analytic  structure  of  the  solution.  The  appli¬ 
cation  is  therefore  rather  blind,  so  the  result 
should  be  viewed  with  skepticism.  How¬ 
ever,  the  pattern  of  convergence  is  often 
manifested  so  convincingly  that  it  speaks 
for  the  accuracy  of  the  tinal  results. 

Consider  the  Cj  transformation,  which  is 
the  simplest  one.  If  three  partial  sums  Vb 
S„.  and  5,^.2  of  a  series  are  known,  then 

^  1  ■5n+ 1  ^n-l  ~ 

e,  - - - — 

5n+l+  ^n-1-  25n 


Figure  28.  Eulerized  series  for  inward  spherical 
solidification. 

(385) 


The  success  of  the  Cj  transformation  in  improving  the  convergence  results  because  it  yields  the  exact 
sum  if  applied  to  a  geometric  series.  Therefore,  it  works  best  on  series  with  nearly  geometric 
coefficients. 

To  illustrate  the  application  of  Shanks  transformation,  we  consider  the  three-term  strained 
coordinate  solution  for  t,  eq  250.  Rewrite  eq  250  as 

X  =  To+  eti  +  £2x2  (386) 

where 


Xo  =  ^(l-i|f^)+  J-\|f2ini|/ 

4  2 

21  nv 

15{1-i|/2)2  +  21  (l  -  y*)lni|/+  12  (l  +  i|r*)ln?i|/  +  3  (l  -  \|/^)ln?i|/ 

96^2  In'*!!! 
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Figure  29.  Shanks-transformed  solution  for 
inward  cylindrical  solidification. 


If  we  note  that  =  Tq,  i'n  =  Tq  +  exi  and  5„+i  =  to  +  eti+  8^t2,  the  application  of  Shanks 
transformation  (eq  38S)  to  eq  386  gives 

e  1  =  t*  (Shanks)  = 

t,-et2 

Figure  29  compares  the  Shanks  transformed  solution  given  by  eq  387  with  the  numerical  results  of 
Tao  (1967).  The  agreement  is  exceptionally  good  even  at  e  =  3.  This  is  a  big  improvement  over  the 
predictions  of  eq  386  which  gives  valid  results  only  up  to  e  =  0.8  (Asfar  et  al.  1979). 


10.  CONCLUDING  REMARKS 

This  review  has  demonstrated  the  usefulness  of  perturbation  techniques  to  analyze  heat  transfer 
problems  involving  freezing  and  melting.  The  perturbation  approach  has  proved  effective  and 
convenient  in  one-dimensional  situations,  and  thus  the  discussion  was  mostly  confined  to  these 
situations.  However,  the  method  has  also  enjoyed  limited  success  with  two-dimensional  cases;  such 
studies  have  been  briefly  mentioned  here  but  furtherdetails  can  be  found  in  the  appropriate  references. 
Despite  their  limited  success  in  more  complex  problems,  perturbation  methods  often  prove  invaluable 
in  illuminating  the  physics  of  the  problem. 

This  monograph  has  been  written  to  serve  two  purposes.  The  first  purpose  is  to  assist  the  unfamiliar 
reader  in  undei'standing  the  perturbation  techniques  and  seeing,  with  the  help  of  detailed  mathematics, 
how  these  techniques  are  applied  to  specific  problems.  The  second  purpose  is  to  bring  together  in  a 
single  publication  the  essence  of  a  large  body  of  literature  on  the  subject  so  that  it  can  serve  as  a 
reference  for  future  studies. 
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